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(4)  INTRODUCTION 

The  “Breast  Cancer  Screening  Using  Photonic  Technology”  research  project  is  devoted  to 
developing  techniques  for  obtaining  direct  two-dimensional  (2-D),  and  tomographic  three- 
dimensional  (3-D)  images  of  cancerous  lesions  of  human  breast  that  make  using  noninvasive 
near-infrared  (NIR)  light.  The  imaging  method  involves  illuminating  the  specimen  with 
ultrashort  NIR  pulses  of  laser  light  and  construction  of  images  using  two  approaches.  The  direct 
or  the  shadowgram  method  utilizes  the  image  bearing  component  of  the  forward-transmitted 
light,  while  the  inverse  reconstruction  method  makes  use  of  the  measured  transmitted, 
forward-scattered  or  backscattered  light  intensity  profiles,  known  or  estimated  optical  properties 
of  the  sample,  a  model  for  light  propagation  through  turbid  media  and  a  sophisticated  computer 
algorithm  to  construct  images  of  the  interior  structure  of  the  specimen. 

Significant  advances  in  developing  both  of  these  approaches  were  made  during  the  current 
reporting  period.  Correlating  the  results  with  that  obtained  from  available  methods  assesses 
efficacy  of  the  methods  so  developed.  We  have  initiated  correlating  the  results  of  optical 
measurements  with  that  obtained  from  nuclear  magnetic  resonance  (NMR)  on  the  same  samples. 

(5)  BODY 

The  tasks  performed  and  the  progress  made  during  the  current  reporting  period  may  be 
broadly  grouped  as  follows: 

5.1  Time-sliced  and  spectroscopic  in  vitro  imaging  of  human  breast  tissues, 

5.2  Correlation  with  NMR  results, 

5.3  Analytical  solutions  of  photon  transport  equation,  and 

5.4  Development  of  the  inverse  image  reconstruction  method. 

We  will  briefly  outline  our  accomplishments  in  each  of  these  areas,  and  refer  to  appended 
publications  (Appendices  1-4)  for  detailed  description  where  applicable. 

5.1  Time-sliced  and  Spectroscopic  in  vitro  Imaging  of  Human  Breast  Tissues 

We  pursued  the  time-sliced  and  spectroscopic  direct  2-D  imaging  of  in  vitro  normal  and 
cancerous  human  breast  tissues  (TO  5-7,  Tasks  12-15)  using  the  experimental  arrangements 
developed  and  improved  during  the  earlier  reporting  periods.  The  time-sliced  imaging 
arrangement  uses  -150  fs,  1  kHz,  800  nm  pulses  from  a  Ti:sapphire  laser  and  amplifier  system 
for  sample  illumination  and  an  ultrafast  gated  intensified  camera  system  (UGICS)  for  recording 
2-D  images  using  different  temporal  slices  of  the  transmitted  light  pulse. 

We  now  routinely  use  UGICS  for  recording  time-sliced  images  after  a  comparison  (TO  3, 
Task  7)  of  the  UGICS  and  optical  Kerr  gate  (OKG)  revealed  that  the  UGICS  is  better  suited  for 
the  task.  OKG  provides  gate  width  limited  by  the  response  time  of  the  Kerr  medium  or  the  pulse 
width  and  could  attain  time  resolution  of  a  few  picoseconds,  compared  to  the  -80  ps  gate  width 
provided  by  the  UGICS.  The  UGICS  is  easier  to  operate,  does  not  need  the  demanding  overlap 
of  the  gating  and  signal  beams  as  required  in  an  OKG,  and  the  80-ps  gate-width  is  reasonable 
compromise  for  the  present  application. 

The  spectroscopic  imaging  arrangement  uses  the  1210-1300  nm  output  of  a  modelocked 
Cr'*‘^:forsterite  laser  for  sample  illumination,  a  Fourier  gate  for  rejecting  multiple  scattered  light, 
and  an  InGaAs  near-infrared  area  camera  for  recording  2-D  images,  as  detailed  in  the  two  earlier 
annual  reports. 

Measurements  were  extended  to  more  breast  tissue  specimens  with  infiltrating  ductal 
carcinoma,  and  infiltrating  lobular  carcinoma  from  patients  of  different  ages.  Samples  were 
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obtained  from  our  collaborators  at  the  Memorial  Sloan  Kettering  Cancer  Center  and  National 
Disease  Research  Interchange  under  an  IRB  approval  at  the  City  College  of  New  York.  ^  ^ 

The  results  of  time-sliced  2-D  imaging  experiments  are  consistent  with  our  earlier  results  ’ 
that  images  recorded  with  light  in  the  earlier  time  slices  highlight  the  cancerous  tissues,  while 
those  obtained  with  light  in  the  later  time  slices  highlight  the  normal  tissues.  Time-sliced 
imaging  can  thus  separate  out  normal  and  cancerous  tissues  in  excised  specimens.  Some  of  these 
results  are  also  presented  in  Appendix  1.  A  sequence  of  these  2-D  images  can  be  used  for  3-D 
image  reconstruction. 

We  are  pursuing  multi-wavelength  differential  imaging  on  normal  and  cancerous  breast 
tissues  (TO  5-7,  Task  15)  to  verify  and  build  statistics  on  our  observation  of  wavelength- 
dependent  difference  in  light  transmission  through  the  cancerous  and  normal  tissues  that  was 
first  reported  in  the  last  Annual  Report.  The  pattern  repeats  in  a  majority  of  the  samples  we 
investigated.  We  are  working  on  defining  a  wavelength-dependent  parameter  whose  values 
would  be  indicative  of  the  state  of  tissue  as  normal  or  cancerous. 

5.2  Correlation  with  NMR  Results 

One  of  the  objectives  of  the  project  was  to  correlate  the  results  of  optical  measurements  with 
those  obtained  from  other  commonly  used  methods  (TO  5-7,  Task  17),  such  as,  pathology,  x- 
ray,  and  NMR.  We  routinely  correlate  our  optical  measurements  on  excised  breast  tissue 
samples  with  pathological  evaluation,  which  is  not  suitable  for  in  vivo  measurements.  NMR  can 
provide  an  in  vivo  assessment.  So  we  have  initiated  measurements  on  the  same  excised  breast 
tissue  specimens  using  both  optical  and  NMR  techniques  to  investigate  how  the  results  correlate. 
NMR  measurements  were  carried  out  at  the  Memorial  Sloan  Kettering  Cancer  Center  by  our 
collaborator  Dr.  Jason  Koutcher.  Figures  1  (a)- 1(f)  present  a  photograph  of  the  sample,  a  NMR 
image,  time-sliced  optical  images  at  delay  times  of  100  ps,  and  300  ps,  as  well  as,  corresponding 
spatial  intensity  distributions.  The  100-ps  image  highlights  the  cancerous  tissue,  which  is  in 
good  qualitative  correlation  with  the  NMR  results  displayed  in  Fig.  1(b).  We  will  pursue  these 
measurements  further  on  other  samples  to  obtain  quantitative  correlation. 

5.3  Analytical  Solutions  of  Photon  Transport  Equation 

Development  of  inverse  image  reconstruction  methods  requires  a  theoretical  model  that 
provides  an  accurate  description  of  photon  transport  through  highly  scattering  turbid  media.  We 
have  built  on  the  earlier  successes^’^  of  our  theoretical  endeavor  (TO  4,  Task  10;  TO  8,  task  18), 
and  extended  those  in  two  very  important  ways.  First,  we  have  developed  an  analytical  solution 
of  the  polarized  photon  transport  equation  in  an  infinite  uniform  medium  using  Cumulant 
Expansion  {Appendix  2)  that  can  take  into  account  the  vector  nature  of  electromagnetie  radiation 
and  can  handle  the  polarization  property  of  light.^  This  is  a  highly  significant  advance  over  our 
Cumulant  Solution  of  Elastic  Scalar  Boltzmann  Transport  Equation  in  an  infinite  uniform 
isotropic  medium,^'*  which  in  itself  was  a  major  advance  over  the  commonly  used  diffusion 
approximation  that  fails  to  adequately  account  for  ballistic  and  snake  photons. 

Second,  our  Cumulant  Solution^  '^  of  the  scalar  equation  has  been  extended  to  semi-infinite 
and  slab  geometries  {Appendix  3)?  This  extension  makes  the  solution  more  suited  for  practical 
application  which  always  involve  finite  geometries. 
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Normal 


Cancer 


Cancer  Normal 


(c)  100  ps  (d)  Pixels 


(e)  300  ps 


Fig.  1(a)  Photograph  of  the  specimen  showing  the  normal  fibrous  and  cancerous  (invasive  ductal  carcinoma)  tissues 
in  the  specimen.  The  pink  dots  represent  orientation  markings  made  on  the  sample  for  histology  following  optical 
and  NMR  measurements,  (b)  A  Ti  weighted  NMR  image  of  the  specimen  showing  the  normal  and  cancerous  parts, 
(c)  A  time-sliced  optical  image  recorded  for  a  delay  time  of  100  ps,  and  the  corresponding  spatial  intensity 
distribution  (d)  integrated  over  the  horizontal  region  denoted  by  the  dashed  box.  Cancerous  region  is  highlighted  in 
this  early-light  image,  (e)  Another  time-sliced  optical  image  recorded  for  a  delay  time  of  300  ps,  and  the 
corresponding  spatial  intensity  distribution  (f)  integrated  over  the  same  horizontal  region  as  in  (c).  Transmitted 
intensity  is  higher  in  the  later-light  image.  The  zero  time  is  taken  to  be  the  time  of  arrival  of  light  through  the  glass 
cell  containing  clear  water  instead  of  the  tissue. 
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5.4  Development  of  Inverse  Image  Reconstruction  Methods 

We  have  completed  the  development  of  a  time-resolved  Fourier  optical  diffuse  tomography 
approach*  for  inverse  image  reconstruction  (TO  4,  Tasks  8, 9;  TO  8,  Task  18)  that  we  presented 
in  the  last  Annual  Report  and  attach  here  as  Appendix  4.  The  approach  allows  use  of  both 
transmitted  and  backscattered  photons  from  the  scattering  medium  to  obtain  3-D  images  of 
embedded  absorbing  and  scattering  inhomogeneities.  It  introduces  the  concept  of  propagation  of 
spatial  Fourier  component  of  the  scattered  wave  field  inside  the  scattering  medium,  and  then 
develops  a  new  optical  diffuse  imaging  methodology  based  on  that  theory.  A  test  run  using 
simulated  data  was  able  to  reconstruct  images  of  four  inhomogeneities  located  up  to  2  cm  below 
the  surface  of  a  human  tissue-like  semi-infinite  scattering  medium  using  backscattered  photons. 
We  are  pursuing  testing  of  this  and  other  methods  developed  in  this  project  for  inverse 
reconstruction  of  objects  using  experimental  data.  An  important  issue  is  a  choice  of  reference 
medium  (TO  4,  Task  10)  and  we  have  used  Intralipid  suspension  of  different  concentration,  and 
suspension  of  Ti02  particles  of  different  size  in  water.  We  are  also  searching  for  matching 
reference  materials  for  breast  tissue. 

We  are  building  forward  model  and  inversion  algorithms  based  on  our  cumulant  solutions  of 
the  scalar  radiative  transport  equation.  Models  using  the  first  two  cumulants  has  been  built  and 
tested  using  simulated  data.  The  model  describes  long-time  behavior  well.  Higher-order 
cumulants  will  be  needed  to  account  for  early-time  behaviors,  and  is  being  actively  pursued. 

(6)  KEY  RESEARCH  ACCOMPLISHMENTS 

•  Verified  earlier  inferences  that  time-sliced  1-D  transillumination  images  recorded 
with  earlier  temporal  slices  of  transmitted  light  highlight  cancerous  tissues  while 
those  recorded  with  later  slices  accentuated  normal  fibrous  tissues  by  extending 
measurements  to  more  breast  tissue  specimens. 

•  Initial  results  show  good  qualitative  correlation  between  results  of  time-sliced  and 
NMR  imaging  methods. 

•  Carried  out  spectroscopic  imaging  experiments  on  more  breast  tissue  specimens  and 
trying  to  verify  if  the  ratio  R  of  light  intensity  transmitted  through  the  cancerous 
tissue  to  that  through  the  corresponding  normal  tissue  as  a  function  of  wavelength 
could  be  used  as  a  useful  parameter  for  cancer  identification. 

•  Developed  analytical  solutions  of  the  polarized  photon  transport  equation  that  enable 
description  of  the  polarization  properties  of  light  transmitting  through  a  highly- 
scattering  medium. 

•  Extended  Cumulant  Solution  of  scalar  transport  equation  to  semi-infinite  and  slab 
geometries  making  it  more  suited  for  practical  application  that  involves  finite 
boundaries. 

•  Developed  a  forward  model  and  computer  algorithm  for  inverse  image  reconstruction 
scheme  that  can  use  both  forward  transmitted  and  backscattered  light  to  provides  fast, 
noise-resistant  3-D  images  of  scattering  and  absorptive  inhomogeneities  at  various 
depths  inside  a  scattering  medium. 

(7)  REPORTABLE  OUTCOMES 
Articles 

1.  S.  K.  Gayen,  M.  Alrubaiee,  M.  E.  Zevallos  and  R.  R.  Alfano,  “Temporally  and  spectrally 
resolved  optical  imaging  of  normal  and  cancerous  human  breast  tissues,”  in  the  Proceedings 
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of  the  Inter-Institute  Workshop  on  In  Vivo  Optical  Imaging  at  the  NIH,  Amir  Gandjbakhche, 
ed.  (Optical  Society  of  America,  Washington,  DC,  2000),  pp.  142-147. 

2.  W.  Cai,  M.  Lax,  R.  R.  Alfano,  “Analytical  solution  of  the  polarized  photon  transport 
equation  in  an  infinite  uniform  medium  using  cumulant  expansion,”  Phys.  Rev.  E  63, 
016606-1(2001). 

3.  M.  Xu,  M.  Lax,  R.  R.  Alfano,  “Time-resolved  Fourier  optical  diffuse  tomography,”  J.  Opt. 
Soc.  Am.  A  18, 1535  (2001). 

4.  A.  Carpenter,  W.  Cai,  M.  Lax,  R.  R.  Alfano,  “Cumulant  approximation  of  the  radiative 
transfer  equation  for  photon  density  in  semi-infinite  and  slab  geometries,”  J.  Opt.  Soc.  Am.  A 
(submitted). 


(8)  CONCLUSIONS 

The  work  carried  out  during  this  reporting  period  builds  on  and  affirms  some  of  our  earlier 
inferences  and  leads  to  some  new  conclusions.  First,  time-sliced  2-D  transillumination  images 
recorded  with  earlier  temporal  slices  of  transmitted  light  highlighted  cancerous  tissues  while 
those  recorded  with  later  slices  accentuated  normal  fibrous  tissues.  Second,  results  of  time-sliced 
measurement  are  in  good  agreement  with  NMR  results  on  same  samples.  Third,  wavelength 
dependence  of  a  ratio,  R  of  light  intensity  transmitted  through  the  cancerous  tissue  shows  the 
potential  to  be  used  as  a  useful  parameter  for  cancer  identification.  Fourth,  analytical  solutions 
of  the  polarized  photon  transport  equation  using  Cumulant  Expansion  that  we  obtained  enables 
description  of  polarization  characteristics  of  light  emerging  from  a  scattering  medium  and  will  be 
useful  for  reconstruction  of  polarization  gated  images.  Fifth,  the  theoretical  formalism  and 
computer  algorithm  for  inverse  image  reconstruction  scheme  using  both  forward  and 
backscattered  light  shows  (with  simulated  data)  the  potential  to  provide  fast  3-D  images  of 
scattering  and  absorbing  objects  at  various  depths  inside  a  scattering  medium. 

We  plan  to  build  on  these  developments  and  pursue  in  vivo  breast  imaging  on  volunteers 
following  testing  of  the  techniques  on  model  breast  structures  made  from  excised  tissues,  and 
have  submitted  a  proposal  to  USAMRMC  seeking  support  for  the  work. 

“So  What  Section” 

The  implication  of  the  first  conclusion  is  that  time-sliced  imaging  offers  the  possibility  of 
highlighting  cancerous  lesions  in  human  breast.  The  second  conclusion  points  to  the 
complementary  nature  of  optical  and  NMR  imaging,  and  the  possibility  of  using  them  together 
for  better  specificity.  The  third  conclusion  indicates  the  diagnostic  potential  of  optical  imaging 
based  upon  multi-spectral  measurements.  X-ray  mammography,  most  often  used  method, 
cannot  diagnose  cancer.  The  fourth  and  fifth  conclusions  together  present  the  possibility  of 
developing  robust  3-D  inverse  image  reconstruction  formalisms,  that  in  addition  to  being 
applicable  for  optical  mammography,  will  be  useful  for  imaging  tumors  in  other  body  organs, 
such  as,  prostate  and  thyroid  glands,  as  well  as,  objects  inside  other  types  of  scattering  media, 
such  as,  cloud,  fog,  smoke,  and  murky  water. 
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Abstract:  Time-sliced  and  spectroscopic  imaging  approaches  were  used  to  obtain  two- 
dimensional  (2-D)  transillumination  images  of  a  composite  in  vitro  human  breast  tissue 
sample  comprising  cancerous  and  normal  fibrous  tissues,  adipose  tissues,  and  a  lymph 
node.  Time-sliced  imaging  approach  used  800-nm,  approximately  130-fs  duration,  1  kHz 
repetition-rate  pulses  from  a  Tirsapphire  laser  system  to  illuminate  the  sample,  and  a 
gated  imaging  system  that  provided  a  variable-position,  ~80  ps-duration  electronic  gate  to 
record  time-sliced  2-D  images.  Images  recorded  with  earlier  temporal  slices 
(approximately,  first  100  ps)  of  the  transmitted  light  highlighted  the  lymph  node  and 
cancerous  tissues,  while  the  later  slices  (later  than  300  ps)  accentuated  the  adipose  and 
normal  tissues.  Spectroscopic  imaging  arrangement  made  use  of  1225  -  1300  nm  light 
from  a  chromium-doped  forsterite  laser  for  sample  illumination,  a  Fourier  space  gate  and 
a  polarization  gate  to  sort  out  a  fraction  of  the  image-bearing  photons,  and  an  InGaAs 
area  camera  for  recording  2-D  images.  Marked  enhancement  of  image  contrast  between 
the  adipose  tissue  and  other  tissues  in  the  specimen  was  observed  when  the  wavelength  of 
imaging  light  was  near  resonant  with  the  1203-nm  optical  absorption  resonance  of  the 
adipose  tissue.  Wavelength-dependent  differences  in  relative  light  transmission  through 
the  normal  and  cancerous  tissues  were  observed. 

OCIS  codes:  (170.3880)  Medical  and  biological  imaging;  (170.6920)  time-resolved  imaging;  (290T050) 
scattering,  turbid  media;  (170.6510)  spectroscopy,  tissue  diagnostics;  (170.3660)  light  propagation  in 
tissues;  (999.999)  Optical  mammography;  (999.999)  near-infrared  absorption  spectroscopy  of  tissues; 

(999.999)  spectroscopic  imaging;  (999.999)  time-sliced  imaging. 

Introduction 

Optical  mammography,  imaging  of  the  interior  structure  of  human  breast  with  light,  is  an  active  area  of 
optical  biomedical  imaging  research.[l-4]  Development  of  optical  breast  imaging  modalities  is  of 
interest  for  several  reasons.  Optical  imaging  methods  are  noninvasive  as  no  ionizing  radiation  is 
involved.  Use  of  different  wavelengths  of  light  has  the  potential  to  provide  diagnostic  information.  In 
contrast  with  x-ray  mammography,  light-based  methods  are  as  apt  to  image  dense  breast  of  a  younger 
patient  as  that  of  an  older  patient.  What  is  even  more  important,  inverse  image  reconstruction  methods 
using  time-resolved  or  frequency-domain  optical  measurements  may  provide  three-dimensional  (3-D) 
tomographic  breast  images. [5-7]  The  ability  to  generate  ultrashort  pulses  and  colors  are  two  major 
attributes  of  light  that  may  be  exploited  to  develop  an  imaging  modality  with  diagnostic  ability. 

In  this  article,  we  present  the  results  of  time-sliced[7]  and  spectroscopic [8]  2-D  transillumination 
imaging  measurements  on  excised  human  female  breast  tissue  specimens  comprising  normal  and 
cancerous  tissues.  Time-sliced  imaging  makes  use  of  different  temporal  slices  of  the  transmitted  light  to 
form  2-D  images  following  the  illumination  of  the  sample  with  ultrashort  near-infrared  (NIR)  pulses  of 
light.  The  thrust  of  the  spectroscopic  imaging  experiment  is  to  examine  if  a  spectroscopic  difference 
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may  be  used  to  enhance  image  contrast,  distinguish  between  different  types  of  tissues  in  a  specimen,  and 
obtain  diagnostic  information. 

Methods  and  Materials 

The  time-sliced  imaging  arrangement  used  800-nm,  approximately  130-fs  duration,  1  kHz  repetition-rate 
pulses  from  a  Ti: sapphire  laser  and  amplifier  system[9]  for  sample  illumination,  and  an  ultrafast  gated 
intensified  camera  system  (UGICS)  for  recording  two-dimensional  images  using  picosecond-duration 
slices  of  light  transmitted  through  the  sample.  The  UGICS  comprised  a  compact  time-gated  image 
intensifier  unit  fiber-optically  coupled  to  a  charge-coupled  device  (CCD)  camera.  It  provided  a 
minimum  gate  width  of  approximately  80  ps  whose  temporal  position  could  be  varied  in  25-ps  steps 
over  a  20-ns  range.  The  average  beam  power  used  in  the  experiment  was  approximately  200  mW.  The 
beam  was  expanded  by  a  beam  expander,  and  a  3-cm  diameter  central  part  of  it  was  selected  out  using 
an  aperture  to  illuminate  the  sample.  The  time-sliced  image  was  recorded  by  the  CCD  camera  and 
displayed  on  a  computer. 

The  experimental  arrangement  for  near-infrared  (NIR)  spectroscopic  imaging  made  use  of  1210- 
1300  nm  continuous- wave  mode-locked  output  of  a  Cr^“‘":forsterite  laser  to  illuminate  the  sample.  A 
Fourier  space  gate[10]  in  conjunction  with  a  polarization  gate[ll]  selected  out  a  fraction  of  the  less- 
scattered  image-bearing  photons  from  the  strong  background  of  the  image-blurring  diffuse  photons.  A 
50  mm  focal-length  camera  lens  placed  on  the  optical  axis  at  a  distance  of  50  mm  from  the  aperture  in 
the  Fourier  gate  collected  and  collimated  the  low-spatial-frequency  light  filtered  by  the  aperture  and 
directed  it  to  the  128x128  pixels  sensing  element  of  an  InGaAs  NIR  area  camera.  The  average  optical 
power  of  the  incident  beam  was  maintained  at  approximately  35  mW  for  all  the  wavelengths  used  in  the 
imaging  experiment  using  appropriate  neutral  density  filters.  The  laser  beam  was  linearly  polarized 
along  the  horizontal  direction. 

The  composite  excised  breast  tissue  sample  used  in  the  experiments  reported  in  this  article  was 
assembled  from  tissues  obtained  following  the  modified  radical  mastectomy  of  a  30-year-old  patient.  It 
comprised  a  lymph  node  (LN)  with  surrounding  tissues,  a  piece  of  adipose  (A)  tissue,  and  a  piece  with 
normal  (N)  and  cancerous  (C)  fibrous  tissue.  Each  of  the  pieces  was  approximately  5  mm  thick,  and 
was  pressed  into  a  5-mm  thick  quartz  cell  to  ensure  uniform  sample  thickness  and  good  optical  contact 
between  the  adjacent  pieces.  According  to  an  accompanying  surgical  pathology  report,  the  cancer  was  a 
poorly  differentiated  carcinoma,  grade  III  with  sarcomatoid  features.  Figure  1(a)  shows  a  photograph  of 
the  exit  face  of  the  composite  sample  wherein  the  locations  of  different  types  of  tissues  are  tentatively 
outlined  and  labeled.  The  tissues  were  made  available  to  us  by  National  Disease  Research  Interchange 
under  an  IRB  approval  from  the  City  College  of  New  York. 

Results 

Time-sliced  Imaging 

Time-sliced  transillumination  images  of  the  sample  for  gate  positions  of  100  ps  and  350  ps  are  displayed 
in  Figs.  1(b)  and  1(c),  respectively.  The  zero  position  was  taken  to  be  the  time  of  arrival  of  the  light 
pulse  through  a  5-mm  thick  quartz  cell  filled  with  water.  The  spatial  intensity  profiles  of  the  images  in 
Fig.  1(b)  and  Fig.  1(c)  integrated  over  two  6-pixel  wide  horizontal  areas  around  the  white  dashed  lines 
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Delay  100  ps 


(b)  LN  A 


Delay  350  ps 


Fig.  l.(a)  A  photograph  of  the  exit  face  (the  side  that  faces  the  camera  in  the  experimental  arrangements)  of  the  composite 
breast  tissue  sample.  LN,  lymph  node;  A,  adipose  tissue;  N,  normal  fibrous  tissue;  C,  cancerous  tissue.  Time-sliced 
transillumination  images  of  the  sample  for  gate  delays  of  (b)  100  ps,  and  (c)  350  ps.  Spatial  profile  of  the  integrate 
intensity  distribution  along  a  horizontal  area  of  6  pixel  vertical  width  around  the  dashed  white  line  that  passes  through  the 
normal  fibrous  tissue  (dashed  line),  or  the  cancerous  tissue  (solid  line)  for  a  gate  delay  of  (d)  100  ps,  and  (e)  350  ps. 

are  presented  in  Figs.  1(d)  and  1(e),  respectively.  The  two  areas  were  chosen  such  that  one  included  the 
normal  fibrous  tissue  in  the  upper  right  part  of  the  sample  while  the  other  included  the  cancerous  tissue 
in  the  lower  right  part  to  enable  close  comparison.  The  time-sliced  100-ps  gated  image  clearly 
highlights  the  lymph  node,  adipose,  normal  fibrous,  and  cancerous  tissue  regions.  The  contrast  is  the 
highest  between  the  lymph  node  that  appears  the  brightest  and  the  adipose  tissue  that  appears  dark  in  the 
image.  The  spatial  intensity  distributions  of  the  100-ps  image,  displayed  in  Fig.  1(d),  show  the  highest 
peak  in  intensity  values  in  the  lymph-node  region  and  a  marked  trough  in  the  adipose  tissue  region 
indicating  much  higher  light  transmission  through  the  lymph  node  and  much  lower  transmission  through 
the  adipose  tissue  region  at  early  time.  More  interesting  is  the  contrast  between  the  cancerous  and 
normal  tissues  in  the  100-ps  image.  As  seen  in  the  right  side  of  the  image  and  the  spatial  intensity 
profiles  of  Fig.  1(d),  light  transmission  through  the  cancerous  tissue  is  significantly  higher  than  that 
through  the  normal  tissue. 

A  markedly  different  situation  is  observed  in  the  350-ps  gated  image  of  Fig.  1(c)  and  the 
corresponding  spatial  intensity  profiles  of  Fig.  1(e).  The  adipose  tissue  region  appears  the  brightest,  and 
the  spatial  intensity  profile  peaks  in  the  adipose  tissue  region  indicating  much  higher  light  transmission 
through  the  adipose  tissue  compared  to  transmission  through  other  tissues  in  the  sample  at  this  later 
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time.  What  is  even  more  noteworthy,  the  difference  in  light  transmission  through  the  normal  and 
cancerous  regions  that  appeared  so  prominent  in  the  profiles  of  Fig.  1(d)  is  not  appreciable  at  this  late 
time.  It  is  reflected  by  the  close  overlapping  of  the  two  profiles  in  the  regions  of  the  normal  and 
cancerous  tissues  in  the  profiles  of  Fig.  1(e).  At  intermediate  times  (not  shown  in  figures)  relative  light 
transmission  through  the  lymph  node  and  cancerous  tissues  decreased  while  that  through  adipose  and 
normal  fibrous  tissues  increased  with  time.  Summarizing  the  time-dependent  transit  of  light,  we  find 
that  the  light  transits  fastest  through  the  lymph  node,  followed  by  that  through  cancerous  fibrous  tissue, 
normal  fibrous  tissue,  and  the  adipose  tissue.  Lower  scattering  or/and  higher  absorption  of  light  by  a 
particular  type  of  tissue  may  make  its  temporal  profile  peak  at  an  earlier  time.  Since  there  is  no  known 
absorption  of  800-nm  light  by  breast  tissues,  we  attribute  the  observed  time-dependent  differences  in  the 
relative  light  transmission  to  the  differences  in  the  light  scattering  characteristics  of  these  tissues.  Our 
results  suggest  that  lymph  node  scatters  the  least,  followed  by  cancerous  tissue,  normal  fibrous  tissue, 
and  the  adipose  tissue. 

These  results  demonstrate  that  time-sliced  imaging  can  highlight  different  types  of  tissues  in  a 
sample.  What  is  even  more  important,  it  can  highlight  normal  fibrous  tissues  from  poorly  differentiated 
carcinoma  (grade  III)  with  sarcomatoid  features.  It  should  be  noted  that  more  pronounced  difference  in 
light  scattering  characteristics  between  normal  fibrous  tissues  and  infiltrating  ductal  carcinoma  was 
observed[12]  than  that  observed  in  this  study  between  normal  fibrous  tissues  and  poorly  differentiated 
carcinoma  (grade  III)  with  sarcomatoid  features.  While  both  types  of  tumors  scatter  less  than  normal 
fibrous  tissue,  ductal  carcinoma  tumor  is  even  less  scattering  than  the  poorly  differentiated  carcinoma 
(grade  III)  with  sarcomatoid  features.  Experiments  with  different  types  of  tumors  at  different  stages  of 
development  are  needed  to  generate  quantitative  data  on  light  scattering  characteristics  of  tumors. 

Spectroscopic  Imaging 

A  spectroscopic  difference  between  different  types  of  tissues  in  a  specimen  is  expected  to  provide  some 
distinguishable  signature  in  a  transillumination  image.  In  order  to  test  if  this  signature  may  be  realized  in 
practice,  we  obtained  images  of  the  sample  with  1225-nm  light  that  is  near-resonant  with  the  adipose 
tissue  optical  absorption  resonance  around  1203  nm,[13]  as  well  as,  with  light  of  wavelengths  away 
from  the  resonance.  Figures  2(a)  and  2(b)  show  a  'near-resonant  image  recorded  with  1225-nm  light, 
and  a  typical  'nonresonant  image'  recorded  with  1300-nm  light,  respectively.  Figures  2(c)  and  2(d) 
display  the  corresponding  spatial  intensity  profiles.  The  solid  line  in  each  of  these  figures  shows  the 
profile  integrated  over  a  6-pixel  wide  area  around  the  long  dashed  line  that  runs  the  entire  length  of  the 
corresponding  image  and  includes  the  cancerous  tissue  region  in  the  lower  right  part  of  the  image.  The 
dashed  curve  superimposed  on  the  solid  curve  shows  the  profile  integrated  over  a  6-pixel  wide  normal 
fibrous  tissue  area  around  the  upper  smaller  dashed  line  in  the  corresponding  image.  The  solid  and  the 
dashed  curves  in  the  right  side  of  the  profiles  thus  enable  comparison  of  light  transport  characteristics 
through  normal  and  cancerous  tissues  in  the  specimen. 

The  salient  features  of  the  spectroscopic  images  and  corresponding  profiles  are:  (a)  the  adipose 
tissues  appear  much  darker  (less  light  transmission)  than  other  tissues  in  the  near-resonant  1225-nm 
image  as  compared  to  that  in  the  off-resonance  1300-nm  image;  (b)  cancerous  tissues  appear  brighter 
(higher  light  transmission)  than  the  normal  tissues  in  both  the  images;  (c)  while  the  overall  light 
transmission  through  the  normal  region  remains  approximately  at  the  same  level,  that  through  the 
cancerous  region  is  significantly  higher  at  1225  nm  than  at  1300  nm;  (d)  transmission  through  the  lymph 
node  exhibits  a  wavelength-dependent  variation  as  well,  being  higher  at  1300  nm  than  at  1225  nm.  The 
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A,  =1225  nm  ^  =1300  nm 


Figure  2.  Spectroscopic  2-D  transillumination  image  of  the  breast  tissue  sample  described  in  the  text  for  light  of  wavelength  (a)  1225  nm, 
and  (b)  1300  nm.  Corresponding  spatial  profiles  are  shown  in  (c)  and  (d),  respectively.  Spatial  profiles  are  the  integrated  intensity 
distribution  along  a  6-pixel  wide  horizontal  area  around  the  white  dashed  line  covering  the  entire  length  of  the  sample  including  the 
cancerous  tissue  region  (solid  line  in  the  profiles)  and  that  around  the  small  dashed  line  denoting  the  normal  tissue  region  (broken  line  in 
the  profiles). 

change  in  image  contrast  was  the  maximum  for  the  adipose  tissue  as  the  wavelength  of  imaging  light 
was  changed  from  1225  to  1300  nm.  Adipose  tissue  region  appeared  as  a  much  deeper  trough  in  the 
spatial  intensity  profile  of  the  1225-nm  image  compared  to  that  for  the  1300-nm  image.  For  a  more 
quantitative  description  of  the  observed  behavior,  we  monitored  the  image  contrast,  C(X)  =  (Ip  -hVih 
+Ia),  where  is  the  optimal  intensity  value  at  wavelength  A  on  the  spatial  profile  of  the  image  at  the 
adipose  tissue  location,  and  IpiX)  is  the  corresponding  intensity  in  the  immediate  fibrous  tissue  region. 
Value  of  contrast  at  1225  nm  is  0.27  and  0.10  at  1300  nm.  As  the  laser  output  was  tuned  away  from 
1225  nm  to  off-resonance  wavelengths,  the  contrast  between  the  adipose  and  fibrous  regions  in  the 
images  decreased  from  that  maximum  value  of  0.27  towards  0.10.  These  results  clearly  demonstrate 
that  an  appreciable  spectroscopic  difference  may  significantly  enhance  the  contrast  between  different 
types  of  breast  tissues  in  a  transillumination  image,  and  is  consistent  with  our  previous  results  with 
adipose  and  fibrous  human  breast  tissues. [8] 

Even  more  promising  and  interesting  is  the  wavelength-dependent  difference  in  light 
transmission  through  the  cancerous  and  normal  tissues.  As  a  measure  of  this  difference  we  may  monitor 
the  ratio,  R  of  light  intensity  transmitted  through  the  cancerous  tissue  to  that  through  the  corresponding 
normal  tissue.  Taking  the  averaged  intensity  values[14]  around  the  middle  of  the  normal  and  cancerous 
tissues  (Pixel  #110  marked  by  vertical  dashed  line  in  figures  2(c)  and  2(d)),  we  obtain  the  value  of  R  to 
be  1.5  for  1225  nm  and  1.2  for  1300  nm,  a  significant  difference.  We  observed  similar  wavelength- 
dependent  variation  in  R  for  ductal  carcinoma  and  normal  breast  tissue  samples  as  well.  More 
measurements  involving  normal  tissues  and  tissues  with  different  types  and  stages  of  cancer  are  needed 
to  examine  if  R  can  be  a  parameter  whose  values  would  be  indicative  of  cancer. 
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In  summary,  the  spectroscopic  and  time-sliced  imaging  approaches  show  tissue  selectivity.  A 
combined  spectroscopic  and  time-sliced  imaging  approach  along  with  inverse  image  reconstruction 
method  has  the  potential  to  provide  more  information  even  than  the  x-ray  techniques. 
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An  analytical  solution  for  time-dependent  polarized  photon  transport  equation  in  an  infinite  uniform  isotro¬ 
pic  medium  is  studied  using  a  circular  representation  of  the  polarized  light  and  expansion  in  the  generalized 
spherical  functions.  We  extend  our  cumulant  approach  for  solving  the  scalar  (unpolarized)  photon  transport 
equation  to  the  vector  (polarized)  case.  As  before,  an  exact  angular  distribution  is  obtained  and  a  cumulant 
expansion  is  derived  for  the  polarized  photon  distribution  function.  By  a  cutoff  at  the  second  cumulant  order, 
a  Gaussian  analytical  approximate  expression  of  the  polarized  photon  spatial  distribution  is  obtained  as  a 
function  of  the  direction  of  light  and  time,  whose  average  center  position  and  half-width  are  always  exact.  The 
central  limit  theorem  claims  that  this  spatial  distribution  qjproaches  accuracy  in  detail  when  the  number  of 
collisions  or  time  becomes  large.  The  analytical  expression  of  cumulants  up  to  an  arbitrary  high  order  is  also 
derived,  which  can  be  used  for  calculating  a  more  acciuate  polarized  photon  distribution  through  a  numerical 
Fourier  transform.  Contrary  to  what  occurs  in  other  approximation  techniques,  truncation  of  the  cumulant 
expansion  at  order  n  is  exact  at  that  order  and  cumulants  up  to  and  including  order  n  remain  unchanged  when 
higher  orders  are  added,  at  least  as  applied  in  our  photon  transport  equation. 

DOI:  10.1103/PhysRevE.63.016606  PACS  number(s):  42.25.Fx,  42.25 Ja,  42,25.Dd,  42.68.Ay 

L  INTRODUCTION 

Study  of  the  polarized  photon  transport  has  lasted  for 
many  years  since  the  polarized  photon  transport  equation 
(PPTE)  was  formulated  by  Gans  [1]  and  by  Chandrasekhar 
[2].  Recently,  polarization  analysis  of  light  migrating  in  a 
multiple-scattering  medium  has  been  applied  to  broad  fields, 
such  as  diagnostics  of  biological  tissues  [3-5],  atmosphere 
monitoring  [6],  and  communications.  One  goal  is  to  develop 
optical  tomography  with  polarization  analysis  to  enhance 
ability  in  image  reconstruction  of  objects  inside  scattering 
media.  Because  of  the  depolarization  effect  in  a  highly  scat¬ 
tering  medium,  scattered  photons  maintaining  polarization 
are  those  near  ballistic  and  snake  like,  which  suffer  less  mul¬ 
tiple  scattering.  Therefore,  a  tomographic  approach  using  po¬ 
larized  photons  vrill  automatically  exclude  multiple-scattered 
diffusive  photons  which  blur  images.  In  order  to  build  a 
proper  forward  model  for  tomography  using  a  polarization 
analysis,  a  theoretical  study  of  the  propagation  of  polarized 
light  in  scattering  media  becomes  practically  important. 

In  polarized  photon  transport,  the  intensity  of  polarized 
light  scattered  fiom  a  scatterer  along  a  certain  direction  is 
determined  by  many  scattering  processes  at  different  scatter¬ 
ing  planes  consisting  in  a  ray  scattered  from  the  scatterer  and 
rays  incident  to  the  scatterer  from  different  directions.  In 
order  to  properly  describe  this  process,  Kuscer  and  Ribaric 
[7]  employed  a  circular  representation  of  the  polarized  com¬ 
ponents  of  light  and  an  expansion  by  generalized  spherical 
functions  [8]  (or  rotation  matrices  in  angular  momentum 
theory  [9,10]).  The  phase  matrix,  hence,  can  be  analytically 
expressed  by  the  angular  parameters  of  the  incident  and  scat¬ 
ter^  rays  in  fixed  coordinates.  Based  on  this  formalism, 

Herman  and  Lenoble  [1 1]  studied  the  asymptotic  behavior  of 


the  polarized  radiation  field  at  great  depths.  Domke  [12]  con¬ 
structed  a  system  of  singular  eigenfunctions  of  the  PPTE,  for 
which  an  integral  equation  and  then  a  recurrence  relation 
were  derived.  However,  to  our  knowledge,  an  explicit  ana¬ 
lytical  expression  of  the  solution  of  the  PPTE  has  not  been 
obtained.  Numerical  methods,  mainly,  Monte  Carlo  simula¬ 
tions,  are  the  main  tools  in  recent  theoretical  investigations 
of  light  polarization  in  multiple-scattering  media  [5,13]. 

In  this  paper  we  derive  an  analytical  solution  of  the  time- 
dependent  PPTE  in  an  infinite  uniform  medium.  Based  on 
our  results,  inverse  image  reconstruction  of  objects  inside  a 
scattering  medium  using  polarized  light  can  be  developed. 
The  4X4  phase  matrix  is  assumed  to  depend  only  on 
the  scattering  angle  0  in  the  scattering  plane:  P(s,So) 
—  P(s-So)  =  P(cos©).  Under  this  assumption  an  arbitrary 
phase  matrix  can  be  handled.  By  use  of  the  circular  repre¬ 
sentation  of  polarized  light  and  an  expansion  in  the  general¬ 
ized  spherical  function  [7,9],  we  extend  our  approach  in 
solving  the  scalar  (unpolarized)  photon  transport  equation 
[14,15]  using  a  cumulant  expansion  to  the  vector  (polarized) 
case.  Terminating  at  second  order,  an  approximate  Gaussian 
polarized  photon  spatial  distribution  is  obtained  for  a  given 
light  direction  s  as  a  function  of  time  U  Our  solution  for  the 
distribution  in  angle  is  exact,  as  are  the  first  and  second 
cumulants  in  space  at  any  angle  and  time,  which  guarantees 
the  correct  central  position  and  the  correct  half-width  of  the 
spatial  distribution.  After  many  scattering  events  have  taken 
place,  the  central  limit  theorem  claims  that  the  spatial  Gauss¬ 
ian  distribution  calculated  will  become  accurate  in  detail, 
since  all  cumulants  higher  than  the  second  approach  small 
values  relative  to  the  appropriate  power  of  the  second  cumu¬ 
lant.  At  early  times,  the  spatial  distribution  is  narrow:  hence, 
a  distribution  function,  the  mean  position  and  half-width  of 
which  are  exact,  may  provide  an  adequate  description  of  a 
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polarized  beam  in  the  presence  of  noise  and  finite  instrument 
resolution.  An  analytical  expression  of  cumulants  up  to  an 
arbitrary  high  order  is  also  derived  Using  these  higher-order 
cumulants,  through  a  numerical  Fourier  transform,  a  more 
accurate  solution  of  the  PPTE  can  be  calculated 

The  paper  is  organized  as  follows.  Section  n  provides  the 
preliminary  formula  of  the  PPTE,  the  circular  representation 
of  polarized  light  and  the  generalized  spherical  function, 
which  have  been  published  in  previous  literature.  In  Sec.  DI, 
an  exact  solution  of  the  angular  distribution  of  the  polarized 
light  is  derived  and  an  expression  in  the  cumulant  expansion 
of  the  polarized  photon  distribution  is  presented  In  Sec.  IV, 
by  terminating  the  cumulant  expansion  at  second  order  a 
Gaussian  approximate  spatial  distribution  as  a  function  of 
light  direction  s  and  time  t  is  obtained,  and  the  exact  expres¬ 
sions  of  the  first  and  second  cumulants  are  derived.  In  Sec. 
V,  an  expression  of  cumulants  up  to  an  arbitrary  high  order 
is  derived.  A  brief  discussion  and  summary  then  follows  in 
Sec.  VI.  In  Appendix  A,  the  expressions  for  coefficients 
in  Eq.  (15)  are  presented  In  Appendix  B,  analytical 
formulas  for  evaluating  integrals  in  Eq.  (42)  are  derived 

EL  CIRCULAR  REPRESENTATION  AND  GENERALIZED 
SPHERICAL  FUNCTIONS 

In  this  section,  we  summarize  the  description  of  polarized 
light  propagation  in  a  scattering  medium  discussed  in  the 
previous  literature. 

Considering  a  light  beam  traveling  along  a  direction  s,  we 
choose  a  reference  plane  through  the  direction  of  propaga¬ 
tion.  Two  complex  components  of  the  electric  field  E,  such 
as  £:,i  =  aiexp(/^),  the  component  parallel  to  the  reference 
plane,  and  =a2  exp(f5j),  the  component  perpendicular  to 
the  reference  plane,  can  be  used  to  describe  a  single  coherent 
beam.  Four  real  components  were  introduced  by  Stokes  [16], 
each  with  the  dimension  of  the  square  of  a  field  or,  more 
precisely,  an  intensity.  The  four  Stokes  parameters  are  col¬ 
lected  into  a  four-element  array  [17].  The 

component  I  is  the  total  intensity: 

+  (la) 

The  component  Q  describes  a  linear  polarization: 

C=<a;>-(<.i>=(|£,P-|£j">.  (lb) 

The  component  U  describes  a  linear  polarization  45°  relative 
to  the  reference  plane: 

f/=<2a,a2Cos5)=<|£(45")|2-l£:(-45°)P),  (Ic) 

where 

£(±45“)=2"''2(£'ii±£:x). 

The  last  component  V  is  the  difference  between  the  intensi¬ 
ties  of  right-  and  left-circularly  polarized  light: 

V=<2aia2sin^)=(|£j?p-|££,P),  (Id) 

where  the  right-  and  left-circular  components  of  field  are 


In  Eq.  (1),  5=  “  ^2 ;  die  angular  brackets  mean  the  aver¬ 

age  over  many  waves  with  independent  phases  in  a  light 
beam. 

To  give  further  physical  meaning  to  the  symbols,  we  note 
that  the  use  of  the  Cauchy-Schwarz  inequality  leads  to  the 
inequality 

(2) 

For  a  coherent  beam,  which  requires  no  averages  in  Eqs. 
1(a)— 1(d),  the  equality  is  automatically  obeyed.  The  opposite 
extreme  case  is  unpolarized  light  for  which 


and  the  total  intensity  I  is  totally  incoherent.  More  generally, 
the  difference  between  the  left-  and  right-hand  sides  of  the 
inequality,  Eq.  (2),  constitutes  the  incoherent  part  of  the  total 
intensity. 

The  kinetic  equation  for  a  polarized  photon  distribution 
function  I(r,s,r)  as  a  function  of  time  t,  position  r,  and  di¬ 
rection  s,  in  an  infinite  uniform  medium,  ftom  a  point  pulse 
polarized  light  source,  I^^^^(r*“ro)^(s— So)^(t-*0),  is 
given  by  [2] 

^I(  r,  s,  r  )/<?t + cs  •  V  ^(  r,  s,  r) + /4  J(  r,s,  r ) 

P(s,s')[I(r,s',t)-I(r,s,r)]ds' 
+I^o>^(r-ro)^(s-so)^(t-0).  (3) 


The  quantities  in  the  Stokes  parameter  (SP)  representation 
will  be  marked  by  adding  a  superindex,  for  example,  I^. 
With  a  rotation  of  the  reference  plane  through  an  angle  a 
^0  (in  the  counterclockwise  direction,  when  looking  in  the 
direction  of  propagation)  around  the  light  propagation  direc¬ 
tion,  I  varies  as  I'=L(Qf)L  In  the  SP  representation,  the 
relation  is  given  by 


■  /'  ■ 

Q' 

U' 

.  v. 

0  0 

cos  2a  sin  2a 
—  sin2a  cos2a 
0  0 


■  /■ 

Q 

u 

.  V. 

Usually,  a  meridian  plane  containing  the  z  axis  and  the  light 
direction  s  is  used  as  the  reference  plane  for  the  description 
of  the  polarization  state  [2,9]  as  shown  in  Fig.  1.  In  Eq.  (3), 
c  is  the  light  speed  in  the  medium,  is  the  scattering  rate 
(per  unit  time),  fXa  is  the  absorption  rate,  and  P(s,s')  is  a 
4X4  phase  matrix.  The  following  form  of  the  4X4  phase 
matrix  [9]  is  used: 


P(s,s' )  =  L(  7r-;^)P(cos  0)L(  -  A'' )»  (5) 


where  0  is  the  angle  between  light  rays  before  and  after 
scattering,  and  the  matrices  L(  — ;i^')  and  L('7r“-^)  are  those 
required  to  rotate  meridian  planes  before  and  after  scattering 
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FIG.  1.  Geometry  of  the  scattering  plane  and  the  reference 
planes  related  to  the  incident  ray,  and  the  scattered  ray. 

The  dark  plane  is  the  scattering  plane.  is  the  angle  be¬ 
tween  the  meridian  plane  (s^)  and  the  scattering  plane.  is  the 
angle  between  the  meridian  plane  (s',z)  and  the  scattering  plane. 


onto  or  from  a  local  scattering  plane,  as  shown  in  Fig.  1.  The 
intrinsic  property  of  scattering  mechanism  is  described  by 
the  4X4  scattering  function  P(cos0),  which  involves 

COS0=S‘S'. 

It  is  convenient  to  use  a  representation  of  the  polarized 
light  in  which  L(a)  is  diagonal,  rather  than  Eq.  (4).  A  cir¬ 
cular  parameter  representation  (CP)  was  first  proposed  by 
Kuscer  and  Ribaric  [7].  Later,  a  more  precise  definition  of 
the  CP,  which  matches  with  the  initial  definition  of  polarized 
light  in  the  SP  representation  by  Chandrasekhar  [2],  was  pre¬ 
sented  by  Hovenier  and  van  der  Mee  [9].  Hereafter  we  use 
the  definition  of  the  CP  in  Ref.  [9],  which  is  given  by 
=  where  /o=(/+V)/2,  I2 

=  (<2  +  /^/)/2,  and  I.2^iQ-iU)l2,  or  with 


0  1/0 
10  0  1 
10  0-1 
0  1-/  0 

In  the  CP,  a  rotation  of  the  reference  plane  through  an  angle 
a  around  the  light  direction  causes  to  be  multiplied  by 
exp(-/ma).  Notice  that  Iq  and  I-q  actually  have  the  same 
rotational  property.  For  the  phase  matrix,  the  transform 
=  TP^^“  ^  is  given  between  two  representations. 

In  the  CP,  it  is  convenient  to  expand  the  phase  matrix  P^ 
using  generalized  spherical  functions  (GSF's),  The  general¬ 
ized  spherical  functions,  which  are  related  to  irreducible  rep¬ 
resentations  of  the  rotation  group  on  three  nonzero  Euler’s 
angles,  are  defined  as  follows  [8]. 

For  /^sup(|Aw|,|rt|)  and  />t==cos  6, 
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with 

This  function  is  directly  related  [9]  to  the  rotation  matrix 
^  angular  momentum  theory  [10]  by  dLiO) 
=(i)"“'"P^„(cos  0).  Some  symmetry  properties  of 
are  =  The  orthogonaUty  re¬ 

lation  for  Pm^P^  siven  by  [8,9] 

(  -  I)'"'*'"  J_^PL,n(P)^m.«(p)<^P=  57:^7  dw  •  (8) 

The  phase  matrix  in  the  CP  can  be  expressed  using  the  gen- 
eralized  spherical  functions  [7,9],  For  notational  simplicity, 
in  the  following  the  quantities  without  a  superindex  are  un¬ 
derstood  to  be  in  the  CP.  Denoting  s=(ac,<^)  and  s' 
=  the  addition  theorem  of  GSF’s  [9]  is  given  by 

(see  Fig.  1) 

exp(im;k')P^.„(cos0)exp(/«;v') 

/ 

5=  —I 

Xexp[-/5(<A-<^')].  (9) 

Using  this  addition  theorem  of  GSF’s,  the  variables 
and  0  in  Eq,  (5)  can  be  eliminated,  and  the  components  of 
the  phase  matrix  in  the  CP  can  be  expressed  using  the  angu¬ 
lar  parameters  of  the  incident  and  scattered  ray  in  fixed  co¬ 
ordinates.  If  we  expand  elements  of  the  CP  phase  matrix  in 
the  scattering  plane,  P;;,„(cos0),  by  GSF’s, 

P„,„(cos0)=^E  pL7’L.«(cos0), 

then,  using  Eq.  (9),  the  4X4  phase  matrix  in  fixed  coordi¬ 
nates  can  be  written  as 

pL  2^  i-iypLjp)pUp') 

Xexp[-/5(<^-<^')],  (10) 

with  indices  m,  n  =  2,0,— 0,— 2  and  /^sup(|m|,|n|). 

The  coefficients  provide  an  intrinsic  description  of 
the  scattering  mechanism.  In  most  useful  cases,  the  coeffi¬ 
cients  have  the  properties  [7,9]  that  (i)  and  p^-m 
are  real,  (ii)  p^„„=p‘„„=pL„-„ ,  and  (iii)  P2o=[P2-o']* 
asterisk  means  complex  conjugate).  Therefore,  for  each  I 
^2,  there  are  six  independent  real  elements  Pii^  Po-o » 
p^_2,  ReC^y,  Im[p2o]-  For  ^  =  0  or  1,  only  p‘oo  and  p{,_o 
are  nonzero.  These  numerical  coefficients  were  calculated 
using  Mie  theory  for  some  examples  by  De  Rooij  and  van 
der  Stap  [18].  These  ,  together  with  and  are 
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panuneters  that  describe  the  nature  of  the  scattering  process 
and  are  treated  as  known  in  our  solution  of  the  transport 
equation. 


2/+1 


fL  (r)  ^  ^ 

'""o'  477  ^ 


2  [sLjiexp(-kiO,  /=1,2,3,4. 


(15) 


in.  DERIVATION 


with  the  eigenvalues  given  by 


Having  the  above  knowledge,  we  analytically  solve  Eq. 
(3)  in  an  infinite  uniform  isotropic  medium.  Using  a  proce¬ 
dure  similar  to  that  discussed  in  Refs.  [14]  and  [15],  we  first 
study  the  dynamics  of  the  photon  distribution  in  the  light 
direction  space  in  the  CP,  F(s,so,t),  which  is  a  vector  of 
four  components,  on  a  spherical  surface  for  s  of  radius  1.  The 
kinetic  equation  for  F(s,So,0  can  be  obtained  by  integrating 
Eq.  (3)  over  whole  space  r.  The  spatial  independence  of  /tj , 
fia,  and  P(s,s')  retains  translation  invariance.  Thus  the  in¬ 
tegral  of  Eq.  (3)  obeys 


^F(s,So,f)/^t+AiaF(s.So»^) 


F(s,So,t)- J  P(s,s')F(s',so,0<^s' 


=I<®>5(s-So)<5(f-0). 


(11) 


Since  the  integral  of  the  gradient  term  over  all  space  van¬ 
ishes,  as  shown  in  Ref.  [14],  if  we  expand  F(s,So,/)  in 
GSF’s,  its  /  components  should  not  be  coupled  to  each  other. 
The  mth  component  of  F(s,So  ,t),  with  the  initial  polarization 
in  unit  «o  1*®  expanded  in  GSF’s  in  the  following 

form: 


Fmno(s>®0»0~2  Fj„„^(r)2  (  l)^^L.j(^)^J,no(Fo) 

Xexp[  — i’s(<^— <^o)]®^P('"F'a^)’  (^^^ 

with  m,  no=2,0,-0,-2,  /^sup(|OT|,|n|).  When  Sp  is  set 
along  the  z  direction  and  the  initial  reference  plane  is  set  as 
the  x-o-z  plane,  Eq.  (12)  specializes  to 


F„„^(s,z,r)  =  S  Fi,„^(r)FL,«o(F)®xp(  -  ino<^)®xp(  -  FaO- 

(13) 

Substituting  Eq.  (12)  [or  Eq.  (13)]  into  Eq.  (11),  using  the 
expression,  Eq.  (10),  of  the  phase  matrix,  and  the  orthogo¬ 
nality  relation  of  GSF’s,  Eq.  (8),  an  analytically  solvable 
equation  for  F^fnn^t)  for  each  /  is  obtained: 

dF‘  {t)ldt='2  (14a) 

u  „  w 


kH(i/2){(n{K,+n^2±n{,_o±n'2_2) 


(n{«-n'22±nLo+nt2)^ 


+  Re(n2o) 
—Im(n2o) 


211/2^ 


(16) 


for  /=  1,2,  and  for  1=3,4,  the  sign-f-before  square  brackets 
in  Eq.  (16)  is  replaced  by  -.  The  constant  coefficients 
[^mnoli  ®an  be  anrdytically  determined  using  standard  linear 
algebra  fiom  the  initial  condition,  Eq.  (14b).  A  detailed  ex¬ 
pression  for  is  presented  in  Appendix  A. 

Equation  (12)  [or  Eq.  (13)],  combined  with  Eqs.  (15)  and 
(16)  and  the  coefficients  [fiLoli  “  Appendix  A,  provides  an 
exact  CP  solution  in  the  light  direction  space.  In  the  SP  rep- 
resentation,  we  have  -  _ 


FSP(s,So,0=T-‘F(s,So,f)T.  (17) 

It  can  be  proved  that  all  components  of  F®‘’(s,So,f)  are  real 
numbers.  The  mth  component  {m= I, Q,U,V\  of  the  angular 
distribution  function  in  the  SP  representation,  with  the  initial 
polarized  state  is  obtained  by 

Ff(s,so,0=[F^'*(s,So.OI®'’(">L.  (18) 


Equation  (17)  serves  as  the  exact  Green’s  function  of  polar¬ 
ized  light  propagation  in  the  light  direction  space.  Since  in  an 
infinite  uniform  medium  this  function  is  independent  of  the 
source  position  tq,  requirements  for  a  Green’s  function  are 
satisfied:  especially,  the  Chapman-Kolmogorov  condition  is 
obeyed:  /ds'  FS‘’(s'',s',r-r')F^‘’(s',s,r'-ro)=FS'’(s",s,7 

—to).  In  fact,  in  an  infinite  uniform  medium,  this  propagator 
determines  all  time  evolution  of  polarized  light,  including  its 
spatial  distribution,  because  displacement  is  an  integration  of 
velocity,  cs(r),  over  time.  The  mth  component  [m 
=7,(2, U,V]  of  the  photon  distribution  function  in  the  SP 
representation,  with  the  source  located  at  ro=0, 

the  initial  direction  Sq  ,  and  the  initial  polarization  is 

given  by 

lf{r,s,t)=I^S^r-cj\t')dt'yis{t)-s)^  ,  (19) 


with  =  The  initial  condition 

F„(s,So,/=0)  =  <5„,„jj5(s-So)  and  the  orthogonality  rela- 
tion,  Eq.  (8),  lead  to 

fLo('=0>=  <5m..o(21+  l)/4^-  (14b) 

The  solution  of  Eq.  (14)  can  be  expanded  in  terms  of  eigen¬ 
states: 


where  (•••)„**  means  the  mth  component  of  the  ensemble 
average  in  the  light  direction  space  in  the  SP  representation. 
The  first  S  function  ensures  that  the  displacement  r-0  is 
given  by  a  path  integral.  The  second  S  function  assures  the 
correct  final  value  of  the  direction.  Equation  (19)  is  a  for¬ 
mally  exact  solution,  but  cannot  be  evaluated  directly.  We 
make  a  Fourier  transform  for  the  first  S  function  in  Eq.  (19), 
then  make  a  cumulant  expansion  [19],  and  obtain 
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/f(r,s,/)=F*%,so,0^  Jrfqexpj  ^q-r+S 


-1  K. 


where  T  denotes  time-ordered  multiplication  [20]  and  F^(s,So,t)  is  given  by  Eq.  (18).  In  Eq.  (20)  the  index  c  denotes  a 
cumulant,  which  is  defined  in  textbooks  of  statistics  [21]  and  statistical  physics  [19].  As  for  an  arbitrary  random  variable  A,  we 
have  (A)<,=(A),  (A\=(A^)-<A)(A),  and  a  general  expression  relating  (A*)  and  (A')^ ,  which  is  given  by 


1  /<A)\'.  1  /(A^)  \'2  1  /(A^V'- 


■S{i-ii-2i2 - ni„ - ), 


Hence,  if  (A'),  /=  1,2, . . .  ,k,  have  been  calculated,  (A')^,  i=  1,2, ...  ,k,  can  be  recursively  obtained  and  conversely  [19]. 
The  kth  moment  (the  term  without  index  c)  wfcdch,  according  to  the  cumulant  expansion  theorem,  is  related  to 
Jdrrj^'"rjJ^(r,s,t).  This  moment  can  be  evaluated  using  a  standard  time-dependent  Green’s  function  approach,  which  is 

given  by 


•  ^  m 

XFSP(sW,s(*“'^rfc-rt-l)sj‘Jl*^•••F*’’(s<^^s(‘^f2-^l)•sjJ^F®’’(s^‘^So>^l-0)I®’’^”^  +(perm.)|,  (22) 

Jm  .  J 


where  the  abbreviation  “perm”  means  all  /:!  —  1  terms  oh- 
tained  by  permutation  of  {//},  /=1, . , .  from  the  first 
term. 

In  Eq.  (22),  is  given  by  Eq.  (17). 

Since  Eq.  (22)  is  obtained  using  a  Green’s  function  approach 
without  making  any  approximation  and  Eq.  (17)  is  an  exact 
expression  of  the  angular  Green’s  function,  Eq.  (22)  provides 
an  exact  formula  for  the  kth  moment.  If  we  are  able  to  ex¬ 
actly  evaluate  Eq.  (22)  up  to  kih  order,  through  Eq.  (21),  we 
can  obtain  the  exact  cumulants  of  the  distribution  up  to  the 
kth  order. 

IV.  GAUSSIAN  APPROXIMATION 
OF  THE  DISTRIBUTION 

Terminating  Eq.  (20)  at  second  order  of  the  cumulant  and 
setting  s  in  Cartesian  coordinates,  integration  over  q  in  Eq. 
(20)  can  be  analytically  performed,  which  leads  to  the  fol¬ 
lowing  Gaussian  approximation  expression  of  the  polarized 
photon  distribution.  When  the  initial  Sq  is  set  along  z,  it  is 
given  by 

/®'’(r,s,r)=  [detZ)^^]''"®^P.“4[(^”^ 


vnthm=l,Q,U,V and  a,  /3=jc,y,z.  In  Eq.  (23),  (Ra)^ rep- 
resents  the  position  of  the  average  center  of  the  distribution, 
and  is  related  to  the  half-width  of  the  spread  of  the 

distribution,  which  is  given  by 


[i>f]a^=[</?a/?^>«-(/?«>M<^^>«]/2.  (24) 

(/?„)*’’  in  Eq.  (23)  and  {RaRfi}m  “  Eq.  (24)  can  be  evalu¬ 
ated  using,  separately,  the  first  order  and  the  second  order  of 
Eq.  (22): 


=  >'t[  [  ds'F^^(s,s'd-t')s'^ 

F“(s,z,/)  Jo  J 


XFS‘’(s',z,F)I®‘’<o> 


X(r,-(F,)f)(r^-(F^)f)  , 
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where  (t.c.)  means  that  the  second  term  is  obtained  by  ex¬ 
changing  the  indices  a  and  ^  in  the  first  term.  As  discussed 
at  end  of  last  section,  Eqs.  (25)  and  (26)  provide  exact  ex¬ 
pressions  for  evaluation  of  the  first  and  the  second  moments. 

In  evaluation  of  Eqs.  (25)  and  (26),  it  is  convenient  to  use 
the  components  of  s  in  a  spherical  harmonic  basis: 

,5Q,.y— 

=,[_2‘^'^sin  0^'^*^cos  ^,  +  2'’*'^ sin  (27) 

and  first  calculate  the  corresponding  quantities  in  the  CP. 
Hence  we  write  Eq.  (25)  as 

y.;T-‘<Rj)TI«']  ,  (28) 

J  Jm 

with  a=x,y,z  and  _/=l,0,— 1,  the  indices  of  the  spherical 
haimonic  basis.  1/  is  a  matrix  for  the  transform  from  a 
spherical  harmonic  basis  to  a  Cartesian  basis,  Sa=UajSj, 
given  by 

■_2-i/2  0  2“*^' 

U=  2-^'^i  0  2-^'^i  .  (29) 

0  10. 

<R>  in  Eq.  (28)  is  defined  in  the  CP  as 

ds'2  F„„(s,s’,t-ns'jF„„^(s',z,n, 

"  "  (30) 

where  F„„(S2,s,  ,t2-^i)  «  the  exact  angular  Green’s  func¬ 
tion  in  the  CP,  Eq.  (12).  Similarly,  Eq.  (26)  is  wntten  as 


where  (Ry^Rj,)  is  defined  in  the  CP  as 

X2  ^ )^j2 
«2 


where  y  1  and  j2  sre  spherical  components,  1»  0,  1. 

In  the  evaluation  of  Eqs.  (30)  and  (32),  a  recurrence  rela¬ 
tion  of  GSP’s  is  used,  which  is  directly  derived  from  angular 
momentum  theory  [10].  Defining  Sj= with  y— 1,0, 
“  1,  we  have 

=  7/2  </,l,m,0lZ+/i,m) 

'  h 

X{l,l/t,±j\l+h,n±J)P‘^^!l^j{cos0), 

j,h=  +  U0-l,  (33) 

with  7+1=+*  and  7o~i» 

Clebsch-Goidan  coefficients  in  angular  momentum  theory 
[10],  given  by 


(2/- 1)2/ 


(/+m)(/-/n  +  l)]‘“ 
2/(/+l) 


(/-m)(/+m)|‘'2 

11/2 

m 

[  (2/-1)/  1 

/(/+i). 

(/  +  m)(/+m+l)]*'^ 

■(/-m)(/+m  +  l)' 

(2/-1)2/  J 

2/(/+l) 

(/+/n)(/+»i+ 1) 
(2/+2)(2/+3) 
(/+m+ l)(/~wi+ 1 ) 

(7+0(27+3) 

(/— m)(Z-m+ 1)1’^ 

(2/  +  2)(2/  +  3)  1 


with  the  row 


index  (from  above)  ;  =  1,0,-  1  and  the  column  index  (from  left)  ft  -  1,0,- 1.  These  recurrence  relations  of  GSF  s 
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were  provided  in  Ref.  [8]  with  some  misprints.  Substituting  Eqs.  (12)  snd  (13)  into  Eq.  (30),  snd  using  Eq.  (33)  und  the 
orthogonaUty  relation  of  GSFs,  Eq.  (8),  integrations  over  ds’  and  dr'  in  Eq.  (30)  can  be  analyticaUy  performed.  When  the 
final  direction  s=(^,<^),  we  have 


{Rj)mn=c'2 


X 


h  4n- 

with  «  =  2,0,— 0,— 2,  A=  +  1,0,— 1  and 


2  2  (r)(/-h,l,n,0|/,/i)(/-h,l,wo,-7|/,«o~y)» 


exp(-\j  *r)-exp(-\(r) 


exp(-/4ar),  i,y=l,2,3,4. 


Similarly,  integrations  in  Eq.  (32)  can  also  be  analytically  performed.  We  have 

(t){l-h2,l/i2fi%n2){l-h2,Uo-Ji^-J2\l^no-Ji-j2} 

fwi^Tl^  »/I  J  */*Q  ' 

'!<{,l—h2—hi,i^iSi\l—h2,n\){l—h2—hi,l;no,—j\\l—h2,nQ—i^, 
with  «i,  n2“2,0,— 0,— 2,  /ii,A2~  +  1A“1.  and 


(35) 


(36) 


(37) 


exp(  — *'r)-exp(-Xjr) 


Up  to  now,  algebraic  analytical  expressions  for  the  first 
cumulant  (the  average  center  of  the  distribution)  and  the  sec¬ 
ond  cumulant  (the  half-width  of  spread)  have  been  derived. 
Equations  (36)  and  (38)  involve  the  related  scattering  param¬ 
eters:  fis  and  [defined  after  Eq.  (14)],  through  in 
Eq.  (16)  and  [B^rnn^i  Appendix  A,  and  the  absorption 

parameter,  .  Thus  they  determine  the  time  evolution  dy¬ 
namics.  The  final  light  direction  s  appears  as  an  argument  of 
the  generalized  spherical  harmonics  in  Eqs.  (35)  and  (37). 
Substituting  Eqs.  (18),  (35),  and  (37)  into  Eqs.  (28),  (31), 
and  then  Eq,  (24),  the  first  and  second  cumulants  in  the  SP 
representation  are  obtained  as  functions  of  s  and  t.  The  dis¬ 
tribution  function  of  polarized  light  is  then  expressed  by  Eq. 

(23) ,  with  Eq.  (28)  for  the  average  center  position  and  Eq. 

(24)  for  the  width  of  the  spread.  Equation  (23)  produces  the 
mth  Stokes  component  of  polarized  light  at  position  r,  with 
light  direction  s,  as  a  function  of  time  /,  initialed  by  ro=0, 
So=z,  and  polarized  state  in  an  infinite  uniform  me¬ 
dium. 

It  is  easy  to  reduce  the  above  solution  to  the  scalar  (un¬ 
polarized)  case  by  considering  only  the  /q  component.  Be¬ 
cause  (/,l,0,0|/,0)=0  in  Eq.  (34),  Eqs.  (35)  and  (37)  can  be 
greatly  simplified.  Also,  Eq.  (15)  is  reduced  to  (2/+l)exp 


exp(-\j  ^20-exp(-\j/) 


>)  ()^-xy^>)(xj-x;.-''^)j 


l~h'} 


(38) 


I -  ■ 

(— Il{)(/)/47r  in  the  scalar  case.  Notice  that  the  associated 
Legendre  function  =  (/)"*[(/  +  m) !/(/ 

“m)l]*^Po^(/A);  our  formula  reduces  to  that  given  in  Ref. 
[14]  in  the  scalar  case. 

V.  DISTRTOUnON  FUNCTION  ACCURATE  UP  TO  AN 
ARBITRARY  HIGH-ORDER  CUMULANT 

In  order  to  calculate  the  polarized  photon  distribution 
function  with  accuracy  up  to  an  arbitrary  high  order,  it  is 
more  convenient  to  set  all  spatial  and  angular  vectors  in  the 
spherical  harmonics  basis,  similar  to  Eq.  (27),  and  to  evalu¬ 
ate  Eq.  (22)  via  the  CP: 

=  , . r;  [T-  •gov-  ,...,ji  , 

(39) 


with  ;i . jk=  1,0,- 1  and  G(Jk  ,-,ji  ,0  given  by 
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XF(s,s<'=>,r-/i)sj*^F(s<*>,s<*-‘>,4-f*_i)jJ;/^-"F(s(2>,s<'>,r2-fi)^jJMs<‘>,So,/,-0)+penn.|, 

(40) 

where  /,_i)  is  given  by  Eq.  (12).  Using  the  GSF  recurrence  relation,  Eq.  (33),  and  the  orthogonality 

relation  of  GSF’s,  Eq.  (8),  the  integrals  overds<*^-  in  Eq.  (40)  can  be  analytically  performed.  We  obtain,  when  the  initial 
Sq  is  along  z  and  the  final  s=(ff,<^),  that 


[G(y„....yi,/)L„„=  E  Pl,„^_jt^^^^(cos5)exp 


-/|  no- 2  Jf]  <f> 


n  Yh 
/=! 


XT'  x?  2(/-2*=iA/)+1 

z  "‘2j  h  -"h - 7Z - 

"I  1‘k  '>1 


j-i  ^ 

I  h/._y+ 1  1 


/=' 


/  «  *-« 

x(  /+I*i»^0  ^E)  J/’  Jk—g+l 

with  n/=2,0,-0,-2  and  fi/=  1,0,- 1,  /=  1,2, . . .  ,fc,  with 


*-l  K-g+l  \  1 

/“Ej  **-/+!. «0“  J/|| 


*-?+! 


j+perm. 


(41) 


x£dr*J\r*_i*"£"dr,exp[-\j^^i(r-rt)]exp[-\|“''*(rt-tfc-i)]-"exp[-\j^^^=‘''*“^"^'(U-0)]- 


(42) 


Note  that  all  ensemble  averages  have  been  performed.  Equa¬ 
tion  (42)  involves  integrals  of  exponential  functions,  which 
can  be  analytically  performed.  An  explicit  expression  for 
evaluating  integrals  in  Eq.  (42)  is  presented  in  the  Appendix 
B.  Equation  (42)  involves  all  related  scattering  and  absorp¬ 
tion  parameters  and  determines  the  time  evolution  dynamics. 
The  final  direction  of  light,  s,  appears  as  an  argument  of 
GSF’s  in  Eq.  (41).  Substituting  Eq.  (42)  into  Eq.  (41), 
through  Eq.  (39),  which  transfers  to  the  SP  representation 
and  introduces  the  initial  polarized  condition,  and  using  a 
standard  cumulant  procedure,  the  cumulants  as  functions  of 
angle  s  and  time  t  up  to  an  arbitrary  A:th  order  in  the  SP 
representation  can  be  recursively  obtained.  The  final  position 
r  appears  in  Eq.  (20),  and  its  components  can  be  expressed 
on  a  spherical  harmonics  basis,  similar  to  Eq.  (27).  Then, 
performing  a  numerical  three-dimensional  inverse  Fourier 
transform  over  q,  an  approximate  distribution  function 
/^^(r,s,t)  in  the  SP  representation,  accurate  up  to  ^th  cumu¬ 
lant,  can  be  calculated. 

VI.  DISCUSSION 

In  Sec.  Ill,  we  derived  an  explicit  expression  of  the  po¬ 
larized  photon  distribution  function,  which  guarantees  the 
exact  average  central  position  (the  first  cumulant)  and  the 
exact  width  of  spread  (the  second  cumulant).  Moreover,  with 


I 

an  increase  of  collision  events  or  time,  the  distribution  ap¬ 
proaches  accuracy  in  detail  since  the  higher  cumulants  be¬ 
come  relatively  small  compared  to  the  appropriate  power  of 
the  second  cumulant.  If  we  examine  the  spatial  displacement 
after  each  collision  event  as  an  independent  random  variable 
Ar,*,  the  total  displacement  is  XAr,-  (i  =  1, . . .  ,A0,  with  N 
the  number  collision  events,  which  can  be  estimated  by 
If  we  define  Y=(A0”^^2Ar,  ,  the  central  limit  theo¬ 
rem  claims  that  if  V  is  a  large  number,  then 
— n^3.  Therefore,  the  sum  of  N  variables  will  have 
an  essentially  Gaussian  distribution.  At  early  times,  the  pho¬ 
ton’s  spread  is  narrow:  hence,  in  many  applications  the  de¬ 
tailed  shape  is  less  important  than  the  correct  position  and 
correct  narrow  width  of  the  beam,  because  of  the  finite  reso¬ 
lution  of  detection  devices.  In  case  a  more  accurate  distribu¬ 
tion  at  early  times  is  needed.  Sec.  IV  provides  formulas  for 
analytically  calculating  the  higher  cumulants  up  to  an  arbi¬ 
trary  itth  order.  Then,  performing  a  numerical  three- 
dimensional  Fourier  transform,  the  distribution  function  ac¬ 
curate  up  to  the  kth  order  cumulant  approximation  can  be 
obtained. 

In  summary,  we  present  an  analytical  solution  of  the  time- 
dependent  polarized  radiative  transport  equation  in  an  infi¬ 
nite  uniform  isotropic  medium.  The  Green’s  function  for  the 
angular  part  is  exact.  Using  a  cumulant  expansion,  we  can 
analyticily  calculate  the  spatial  cumulants  up  to  an  aibitraiy 
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high  order.  By  terminating  at  the  second  order,  we  have  de¬ 
rived  an  explicit  expression  of  the  polarized  light  distribution 
fiincdon.  This  expression  is  quantitatively  accurate  up  to  the 
second  order  cumulant  approximation.  Namely,  the  center 
position  and  the  half-width  are  always  exact  and  not  modi¬ 
fied  when  higher-order  cumulants  are  added.  The  central 
limit  theorem  claims  that  after  enough  collision  events,  all 
cumulants  higher  than  second  approach  small  values,  and  the 
Gaussian  spatial  distribution  calculated  approaches  accuracy 
in  detail.  Our  results  are  given  in  terms  of  a  distribution  with 
coefficients  that  can  be  calculated  algebraically,  with  moder¬ 
ate  effort  at  the  second  cumulant  level  and  additional  effort 
to  induce  the  third-  and  higher-order  cumulants.  This  analyti¬ 
cal  solution  provides  a  background  distribution  function  for 
further  study  of  optical  tomography  using  polarized  light. 
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the  unique  solution  of  [5m«Q]i  ^  obtained.  For 

given  «o  components  of  construct  a  column 

vector  in  the  space  of  the  direct  product  of  /  X  m .  Combining 
Eqs.  (Al)  and  (A2)  in  iXm  space,  we  obtain  the  following 
matrix  equation: 

AB=C.  (A3) 

A  is  a  16X 16  matrix: 

B  and  C  are  16X1  column  vectors:  Bjxn=[KnQ\j 
C,xm=  Here  A  and  C  are  given,  while  B  is  unknown. 
Equation  (A3)  is  a  standard  form  of  a  group  of  16  linear 
equations.  The  solution  is  given  by 

^iXm“  A/xm^^Ct(A),  (A5) 


APPENDIX  A 

In  this  appendix  we  calculate  ^  (1^)-  Substi¬ 

tuting  Eq.  (15)  into  Eq.  (14),  we  obtain  a  set  of  linear  homo¬ 
geneous  equations 

n 

where  eigenvalues  Xj  (i=l,2,3,4)  are  given  by  Eq.  (16). 
These  equations,  however,  are  not  linearly  independent. 
Adding  the  initial  condition,  Eq.  (14b),  given  by 


with  A;xw  is  obtained  by  replacing  the  (iXm)th  column  in 
the  determinant  of  A  by  the  column  vector  C. 


APPENDIX  B 

In  this  appendix,  we  derive  an  analytical  expression  for 
Eq.  (42)  to  itth  order.  By  defining 

1,  .  .  .  ,k, 

(Bl) 


Eq.  (42)  can  be  written  as 


4  4 


with 


g(bi+b2)t  gbit  j 


Jo  Jo  Jo 


"  bl(F,+b2)(*l  +  ^2+^3)'~  biM^2  +  ^3) 


It  is  easy  to  directly  calculate  Eq.  (B3)  for  a  few  low  k 
orders: 


4*3' 


1 


+ 


(bi+i>2)^2*3  (*l+fr2  +  ^3)(^2  +  ^3)^3 

(B4c) 


I 


(B4a) 


In  each  step  of  integration,  the  difficulty  is  in  determining  the 
constant  term.  In  the  following  we  prove  that  this  term  is 
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given  by 

Equation  (B3)  can  be  written  as  ^  \  ) 


Jo 

(B6) 


(B5) 

Jo 

Using  integration  by  parts  to  Eq.  (B5),  we  obtain  Recursively  applying  Eq.  (B6),  we  obtain 


J 


=  -r- ‘ ^(0 -  .  ,  ,, — 7 

''  ’  bk  bkibk+bk-i) 


+•••+(—  1 
+•••+(-1)*-* 


bk{bu+bk-i)'-'{bk+bk-i^-"+bk-f) 
bk{bk+b,,-yy--{bk+bk-\  +  "-+bi) 


(B7) 


Equation  (B7)  provides  a  formula  to  recursively  evaluate  Eq.  (40)  up  to  fcth  order.  Also,  Eq.  (B7)  produces  the  above- 
mentioned  constant  term.  An  explicit  expression  of  Eq.  (42)  can  then  be  wntten  as 

,  ^  (-l)*exp[2)3,_^+,f] 

Xexp(-X,.^^,t)^^^  ’ 


(B8) 


with  and 


g  J 

j^g,  or  'Z  bf,  j>g, 

J  f—j  fs=:<r-4-l 


J 

2 

/=*+! 


(B9) 


where  bg  is  defined  in  Eq.  (Bl). 
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Cumulant  Approximation  of  the  Radiative  Transfer  Equation  for  Photon 
Density  in  Semi  -  infinite  and  Slab  Geometries 


A.  V.  Carpenter,  W.  Cai,  M.  Lax,  and  R.  R.  Alfano 

Institute  for  Ultrafast  Spectroscopy  and  Lasers 
New  York  State  Center  of  Advanced  Technology  for  Ultrafast  Photonic  Materials  and 

Applications 

and 

Departments  of  Electrical  Engineering  and  Physics 
The  City  College  of  New  York 
and 

The  Graduate  Center  of  the  City  University  of  New  York 
New  York,  NY  10031 


Abstract 

The  cumulant  analytical  solution  of  the  Boltzmann  radiative  transfer  is  extended  for  semi¬ 
infinite  and  slab  geometries  to  provide  a  more  accurate  picture  for  the  transition  time  from  the 
ballistic  to  the  diffusive  regime.  For  the  scattering  structures,  the  method  of  adding  image 
sources  used  in  the  conventional  diffusion  approximation  for  boundary-based  media  is  modified 
and  applied  to  the  cumulant  transport  model.  Numerical  comparisons  between  the  cumulant 
model  and  the  standard  and  the  center-moved  diffusion  models  are  presented.  The  cumulant 
model,  which  reduces  to  the  center-moved  diffusion  model  at  later  times,  gives  a  better 
prediction  of  photon  migration  at  early  times  than  the  other  models.  This  work  is  useful  for 
photon  migration  in  several  applications:  human  tissues,  clouds,  fog,  and  seawater.  Our 
calculations  were  specifically  aimed  at  clouds  and  optical  wireless  communications. 

OCIS  codes:  (030.5620)  radiative  transfer,  (290.1990)  diffusion,  (290.7050)  turbid  media,  and 
(170.5280)  photon  migration,  (010.1310)  atmospspheric  scattering. 

Submitted  to:  Journal  of  the  Optical  Society  of  America  A 
IUSL#-2001  -05 
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I.  Introduction 

Understanding  time-dependent  photon  migration  is  important  for  numerous  light-based 
applications  in  medical,  military,  space,  and  atmospheric  fields.*"^  When  a  light  pulse  enters  a 
scattering  medium,  it  divides  into  ballistic,  snake-like,  and  diffusive  components.’’  The  ballistic 
component  describes  photons  that  travel  virtually  undistributed  through  the  scattering  medium. 
The  snake-like  photons  endure  only  a  few  scattering  events  and  have  only  slight  deviations  from 
a  straight-line  path  through  the  media.  The  diffuse  component,  which  dominants  at  later  times, 
represents  photons  that  have  undergone  multiple  scattering  events. 

The  Boltzmann  radiative  transfer  equation  describes  the  photon  migration  through  a 
scattering  medium.  Due  to  its  difficulty  in  obtaining  a  solution,  scientists  utilize  approximations; 
in  particular,  the  standard  diffusion  approximation  is  used  to  describe  the  transport  of  photons. 
This  approximation  is  based  on  the  assumption  that  photons  are  isotropically  diffuse  from  a  fixed 
source  with  a  constant  diffusion  coefficient  in  a  uniform  medium.  This  is  not  the  case  as  shown 
by  Zevallos.'®  The  diffusion  approximation  is  valid  only  when  the  propagation  distance  is  much 
greater  than  the  transport  mean  free  path,  /, ,  such  as  L  >  6/, ,  where  L  is  the  distance  from  the 
source  to  the  detector.^  The  diffusion  approximation  cannot  accurately  determine  the  distribution 
of  early-arriving  photons  or  when  the  source-detector  distance  is  small.  One  suggested  model 
for  bypassing  the  distance  issue  is  to  allow  the  photons  to  be  diffused  from  an  initial  point 
located  one  transport  mean  length  from  the  source  position  along  the  incident  direction.^  This  is 
referred  to  as  the  center-moved  diffusion  approximation.  While  the  center-moved  approximation 
results  are  better  than  the  standard  approximation  at  later  times,  the  migration  of  the  early 
photons  cannot  correctly  be  described  because  photon  movement  before  1/,  is  neglected.  Other 
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models,  such  as  the  non-Euclidean  diffusion  model,  have  been  introduced  to  bridge  the  gap 
between  the  nondiffusive  and  diffusive  regime  with  limited  success^ 

Most  recently,  cumulant  analytical  expressions  for  the  photon  distribution  function  and 
the  photon  density  as  solutions  to  the  Boltzmann  radiative  transfer  equation  have  been  derived  by 
Cai  and  co-workers  for  an  infinite  uniform  medium.*’  ^  In  the  infinite  environment,  the  cumulant 
approximation  displays  an  accurate  picture  of  the  photon  propagation  in  the  early  time  or 
nondiffusive  regimes.  For  most  practical  applications,  the  photons  propagate  through  a  limited 
region  as  oppose  to  an  infinite  one.  The  effects  of  a  boundary  need  to  be  considered. 

In  this  paper,  the  cumulant  analytical  expressions  are  extended  to  describe  photon 
migration  in  structures  of  semi-infinite  and  slab  geometries.  We  focus  on  photons  traveling  in 
clouds  with  /,  =100m.  Previously,  the  standard  and  center-moved  diffusion  approximations 

have  been  used  to  study  photon  migration  in  the  semi-infinite  and  slab  media  with  various 
boundary  conditions. ”  The  results  from  the  cumulant  transport  approximation  were  compared 
with  the  outcome  from  the  diffusion  approximations.  In  contrast  to  the  diffusion  approximation, 
the  cumulant  expressions  give  a  more  accurate  picture  of  the  photon  migration,  including  the 
ballistic  and  snake-like  phases.  The  boundary  effects  are  viewed  through  amplitude  and  decay 
pattern  variations  in  the  time-resolved  profile  as  the  distance  from  the  source  increases.  At 
longer  times,  the  cumulant  results  approach  those  of  the  center-moved  diffusion  approximation. 
Our  work  will  be  important  for  imaging  and  propagation  of  optical  coded  information  through 
clouds  in  optical  wireless  systems. 
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IL  Photon  density  in  an  inflnite,  uniform  medium 

This  section  describes  the  cumulant  solution  for  an  infinite  uniform  medium.  The  photon 
distribution  function,  l{r,s,t),  or  the  photon  density  function,  Ar(F, /),  which  is  defined  as 

N{r,  t)  =  jds  l{r,s,  t),  can  be  used  to  describe  the  temporal  and  spatial  profiles  of  scattered 


photons  inside  a  turbid  medium.  For  an  infinite  uniform  medium  with  a  short  pulse  originating 
from  a  point  source  located  at  r  =  0 ,  the  temporal  profile  of  the  photon  density  at  position  r 
from  the  source  for  the  standard  diffusion  model  (sDM)  and  the  center-moved  diffusion  model 
(cmDM)  is  given  by 


«(?,()=  W,(r,t)  = 


{ATlDctf' 


exp 


ADct 


ex 


(1) 


with  /„  =  0  for  sDM  and  /„  =  I,  for  cmDM,  where,  s„  is  the  unit  vector  for  the  incident 

m  m  t  ’  ’  o 

direction,  D  is  the  diffusion  coefficient  and  is  given  by  <D  =  l/3[(l - g  is  the 
anisotropy  factor  or  mean  cosine  of  the  scattering  angle,  is  the  absorption  rate  (nsec'‘  ),/i,  is 

the  scattering  rate  (nsec“' ),  and  c  is  the  speed  of  light  (m/nsec ). 

The  photon  distribution  function  and  the  photon  density  in  an  infinite,  uniform  medium 
have  been  derived  by  Cai  et  al  using  an  analytical  solution  of  the  radiative  transport  equation 
using  cumulant  expansion.  ^  A  Gaussian  distribution  with  the  exact  first-order  cumulant,  which 
provides  the  average  center  of  the  distribution,  and  the  exact  second-order  cumulant,  which 
provides  the  time-dependent  diffusion  coefficients,  presents  a  more  accurate  picture  of  the 
photon  migration,  especially  at  early  times.'®  Photon  migration  is  viewed  as  a  photon  cloud 
spreading  anisotropically  from  a  moving  center  with  a  time-dependent  diffusion  coefficient. 
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Based  on  the  cumulant  expansion,''*'  the  photon  density,  N{r,t),  in  an  infinite,  uniform  medium 
with  the  source  located  at  =  0  and  an  incident  direction  along  =  z  is  given  by 


N(r,t)  = 


exp 


iz-R,r 

AD^ct 


exp 


expC-iU^O  (2) 


where,  the  moving  center  of  the  photon,  ,  is  given  by 


8i 


(3) 


The  time-dependent  diffusion  coefficients  are 


D, 


=  :^|-  +  ^f^[l-exp(-w)]+ 

gi  {81-82) 


—r^ - X  [1  -  exp(-  g^r)]  -  [l  -  exp(-  g,t)f 

82{8i-82)  2g, 


D  =D  = 

XX  yy 


:^ |—  +  2/' - X &  -  ®xp(-  «i0]  +  -r-^ - X [1  -  exp(- gjOl  I 

(^l-.?2)  82{8x-82)  J 


(4) 


(5) 


These  equations  are  independent  of  absorption.  They  show  time  is  real  and  time  to  build  up  the 
diffusion  coefficient. 

The  coefficient  g^  is  defined  as 


8i=l^s 


(21 + 1) 


(6) 


where,  is  the  Legendre  coefficient  of  the  phase  function,  P(cos0),  which  is  given  by 

P(cos0)=  (l/A^)^a,P,  (cos$)  and  cos9  =  s  -s^.  Two  special  values  for  g,  are  gg  =  0,  which 
/ 


A3-6 


follows  from  the  normalization  of  the  phase  function  and  gj  =c//,  ,  where  the  transport  mean 

free  path  is  defined  as  =  c/k(i  -  cos  0)1  with  anisotropy  factor  cos  9 . 

Our  analysis  focuses  on  clouds  as  the  test  media.  Figure  1  shows  the  normalized  moving 
center and  normalized  diffusion  coefficients  {D^ll,  mdD^jl, )  for  a  cloud  of  water 
droplets  with  a  refractive  index  of  1.33  uniformly  distributed  in  air  with  the  droplet’s  radius 
equal  to  the  incident  light’s  wavelength.  From  Mie’s  theory  for  light  scattering,  g,  and 
~  0.16;i^  and  ~  0.27  ,  respectively. 

Near  time  /  =  0 ,  the  diffusion  coefficients  (d^,  )  are  approximately  equal  to  zero  and 

the  moving  center  is  moving  at  the  speed  of  light.  This  corresponds  to  the  ballistic  stage  for  the 
photon.  As  time  increases,  the  moving  center  slows  and  the  diffusion  coefficient  increases  from 
zero.  This  region  represents  the  snake-like  phase  of  photon  migration.  In  the  diffusive  regime, 
the  moving  center  stops  at  ll,,  while  the  diffusion  coefficients  approach  lj3 . 

The  two  diffusion  models  have  the  fixed  values  for  Rjl,  at  all  times  (Rjl,  =  0  for  the 
sDM  and  R^/l,  =1  for  the  cmDM)  and  a  fixed  diffusion  coefficient,  D  =  lj3,  as  indicated  in 

Figure  1.  The  results  from  our  model  approach  the  cmDM  results  at  later  times,  while  giving  a 
more  precise  view  at  early  times. 

Based  on  the  cumulant  analytical  expressions,  we  developed  the  cumulant  transport 
model  (cTM)  for  structures.  We  first  examine  the  difference  between  the  cTM  and  the  diffusion 
models  by  computing  the  temporal  profiles  for  the  photon  density  for  each  model  using  the 
infinite  geometry  as  a  point  of  reference.  For  all  three  models,  the  scattering  parameters  are  the 


same  as  that  for  Figure  1. 
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As  displayed  in  Figure  2a,  the  photon  density  peak  for  an  infinite  medium  occurs  sooner 
for  both  the  sDM  and  cmDM  than  for  the  cTM  with  a  source-detector  separation,  =  2/, .  For 

the  cmDM,  the  peak  arrives  much  faster  than  the  speed  of  light.  This  occurs  because  of  the  initial 
placement  of  the  source  one  inside  the  medium.  The  peak  from  the  cTM  is  much  stronger  than 

that  of  the  sDM  due  to  the  exact  cumulants  used  in  our  model.  Decreasing  below  21,  yields 

unrealistic  peaks  for  both  the  sDM  and  cmDM.  For  the  cTM,  the  long  tum-on  time  for  the 
photon  density  is  consistent  with  the  concept  that  no  light  arrives  before  the  ballistic  time. 
As  z,,!^  is  increased,  as  shown  in  Figure  2b  with  z,,,^  ^he  peaks  from  the  cTM  and  cmDM 

become  similar  in  intensity  with  the  cTM  still  displaying  a  longer  tum-on  time  than  the  cmDM. 
As  shown  in  Figure  2c  for  =10/,,  once  we  are  outside  of  the  limits  for  diffusion 

approximation  >  7  /, ),  no  discernible  differences  exist  between  the  cTM  and  the  cmDM. 

III.  Photon  density  in  uniform,  semi-inHnite  and  slab  media 

In  this  section,  we  extend  the  cumulant  solution  for  an  infinite  uniform  medium  to  that  in 
semi-infinite  and  slab  uniform  media  with  vacuum  (or  completely  absorbing)  boundaries.  A 
previously  used  approximate  approach  for  extending  to  a  semi-infinite  medium  is  to  place  an 
image  source  above  the  surface  to  offset  the  real  source.**^’  This  image  source  creates  a 
boundary  condition  A^(r,r)|^  =0at  the  extrapolated  boundary  B  :  at  a  distance,  z^  outside  the 

medium.  The  approximate  value  of  z^  is  given  by'® 


V 


1-r. 


(7) 
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where,  is  the  effective  reflection  coefficient  at  the  physical  boundary.  In  our  case  for 
completely  absorbing  boundaries,  r^=0.  The  value  of  reduces  to  2/3/,  (=0.7/,).  The 
extrapolation  length  can  be  adjusted  to  compensate  for  a  mismatch  in  the  index  of  refraction  at 
the  surface. 

The  photon  density  in  a  semi-infinite  medium  is  given  by 

(8) 

The  form  of  the  photon  density  for  the  image  term,  N ,  (r,r),  is  the  same  as  Ng{r,t),  except  the 
location  of  the  real  source  is  replaced  by  the  image  source  position. 

In  Figure  3a,  an  example  of  the  fixed  locations  for  the  real  source  (r^ )  located  at  (0, 0, 0) 

and  image  sources  ( /; )  for  the  sDM  and  cmDM  are  given.  In  the  case  of  the  cTM,  an  image 
source  is  placed  at  z,  =  -  2(l  +  a)/, ,  as  shown  in  Figure  3b.  At  early  times,  the  distance  from  a 
point  at  the  boundary  B  :  to  the  center  of  N^(r,t)  is  not  equal  to  the  center  of  A/,  (r,r),  so  the 
boundary  condition  iV(r,t)|g  =  0  at  the  extrapolated  surface  is  not  satisfied.  This  is  not  critical 

at  early  times  because  the  photons  move  ballistically  along  in  the  forward  direction  and  the 
diffusion  coefficients  are  near  zero.  Hence,  the  photon  density  at  the  boundary  and  the  boundary 
effect  are  negligible.  After  a  short  period,  t  >  2/,/c,  the  center  of  N,j(r,t)  has  relocated  from 

z  =  0  to  z  =  l,,as  shown  in  Figure  1,  and  the  center  of  N,  (r,t)  moves  to  z,  =  -  (l  +  OfX  •  Then, 

the  extrapolated  boundary  condition  is  satisfied. 

The  extrapolated  boundary  condition  for  a  slab  of  thickness,  L ,  assumed  that  N(z,f)  =  0 
at  two  extrapolated  surfaces  outside  of  the  slab.*'’  This  condition  cannot  be  met  with  the 
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addition  of  a  single  image  source.  To  satisfy  the  extrapolated  boundary  conditions  iV(r,r)|^  =0 
and  A(r,r)|g  =  0  with  Bj :  z  =  -al,  and  B2 :  z  =  L  +  «/,,  an  infinite  series  of  images  dipole 
sources  are  added.  The  solution  now  becomes 


W(r,<)=Ar,(r,l)+  £  S 


(9) 


where,  and  represent  the  photon  density  for  the  positive  and  negative  image 


terms,  respectively.  A^^!r^(r,r)  and  N^'”\r,t)  are  similar  to  N,j{r,t),  except  the  location  of  the 
real  source  is  replaced  by  either  a  positive  or  negative  image  source  located  at  (o,  0,  z" ).  The 

number  of  dipole  sources  that  are  necessary  for  an  accurate  calculation  depends  on  the  optical 
properties  of  the  slab  and  the  maximum  allowed  time  for  the  measurement. 

For  the  sDM  and  cmDM,  the  fixed  locations  for  the  image  dipoles  with  m  =  0  and  m  =  1 
are  given  in  Figure  4a.  Using  the  sDM,  the  locations  of  the  positive  and  negative  image  sources 
for  m  <  0  and  m  >  1  are  z,”  =  2m{L  +  lal,  )for  positive  sources  and  z”  =  2m{L  +  2al,)-  2al, 
for  negative  sources.  For  the  cmDM,  the  locations  for  infinite  dipole  sources  are  similar  to  those 
for  the  sDM.  However,  to  determine  the  source  position  for  the  cmDM,  one  I,  must  be  added  to 

the  positive  source  location  of  the  sDM,  while  one  I,  must  be  deducted  from  its  negative  source 

f  00  00 


location.  For  the  cTM,  in  order  compensate  by  ^ 


the 


locations  of  the  positive  and  negative  image  sources  now  depend  on  whether  they  are  located 
above  or  below  the  slab,  as  shown  in  Figure  4b.  As  with  the  semi-infinite  geometry,  the  centers 
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of  N^{r,t)  and  N^"!l{r,t)  for  the  slab  do  not  offset  each  other  during  the  early  time.  Thus, 
neither  of  the  extrapolated  boundary  conditions,  A^(r,r)|^  =0  orA(r,t)|^^  =0  are  initially 
satisfied. 

The  initial  positions  are  z”  =  2w(L  +  2a/,  )for  positive  sources  and 
z,”  =  2m{L  +  2a/, )  -  2/,  (a  + 1)  for  negative  sources  located  above  the  slab  ( w  <  0 ).  While  for 
dipole  pairs  located  below  the  upper  physical  surface  of  the  slab  ( m  >  0 ),  the  starting  locations 
are  z,™  =  2m(L  +  2a/,)+ 2/,for  positive  sources  and  z,”  =  2m(L  +  2a/, )  -  2/,a  for  negative 
sources.  In  both  cases,  the  movement  of  the  image  dipole  pairs  during  the  early  time  is  towards 
the  two  extrapolated  boundaries  B,  ;z  =  -a/,  and  B2;z  =  T+a/, .  After  a  short  period, 

DO 

t  >  2/,/c,  the  centers  of  N^{r,t)  and  ^A/j'”^(r,r)are  the  same  are  those  of  the  cmDM  for  all 

m  =  -«» 

m .  The  extrapolated  boundary  conditions  are  now  achieved. 

The  results  for  a  uniform  semi-infinite  and  the  slab  media  for  Zj^,,  =  7.0/,  are  shown  in 

Figures  5a  and  5b,  respectively.  In  all  the  slab  geometries,  the  thickness,  L ,  of  each  slab  is  equal 
to  .  For  the  semi-infinite  and  slab  media,  the  early  time  results  are  the  same  as  those  for  the 

infinite  medium  because  the  effects  from  the  flat  surfaces  are  negligible  for  ^  2/, .  For  small 

source-detector  separations  within  either  geometries,  the  cmDM  again  fails  to  yield  a  real  peak 
and  the  peak  from  the  sDM  is  clearly  too  weak  for  any  reasonable  experimental  measurements. 
For  small  z^do » physical  boundary  effects  can  only  be  viewed  far  into  the  diffusive  regime.  As 

z^d^  increases,  the  cTM  and  the  cmDM  yield  nearly  the  same  amplitude.  Because  of  the 
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anisotropic  effects  that  are  retained  by  the  cTM,  the  photon  density  from  the  cTM  is  time-delayed 
from  that  of  the  cmDM. 

The  photon  density  generated  by  the  cTM  for  all  three  geometries  can  be  viewed  for  the 
nondiffusive  =  2.0/, )  and  diffusive  =  7.0/, )  regimes  in  Figures  6a  and  6b,  respectively. 

The  non-diffusive  regime  is  the  same  for  each  geometric  configuration.  More  importantly,  the 
cTM  shows  a  long  tum-on  time  and  does  not  give  a  fake  peak  even  when  is  small.  The 

general  pulse  shape  is  the  same  for  each  configuration,  however,  the  exponential  decay  does  vary 
depending  on  the  media.  The  slab  exhibits  a  faster  exponential  decay,  while  both  the  semi¬ 
infinite  and  slab  media  have  broader  peak  with  a  slower  decay.  As  z,^^  is  increased,  the  damping 

effect  from  reflections  at  the  surfaces  is  evident.  From  a  semi-infinite  medium  to  a  slab  medium, 
the  effective  diffusion  coefficient  is  reduced  as  the  number  of  physical  boundaries  increases.'^ 
Relevant  non-diffusive  information  about  photon  migration  is  present  during  the  early  times  for 
all  three  geometric  configurations.  Our  model,  unlike  the  other  diffusion  models,  does  not 
neglect  or  distort  the  non-diffusive  regime  for  either  the  semi-infinite  or  slab  media  even  when 
the  z,a^  is  small. 

We  used  the  extrapolated  boundary  condition  with  a  non-reflecting  physical 
boundary (r^  =o)  at  z  =  0.  A  cloud  layer  satisfied  the  physical  boundary  condition  with  water 

droplets  located  inside  air  as  the  scatters.  Our  solution  can  easily  be  modified  to  include 
reflecting  boundaries  (r^  ^  o)  by  changing  the  location  of  the  extrapolated  boundaries.  Using  a 

modified  a  created  by  multiplying  the  original  a  (for  =  0)  by  a  factor  (l  +  r^)/(l  -  r^),  we 
can  determine  the  new  position  of  each  image  source.  Paasschens  and  ‘t  Hooft  have  another 
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formula  for  reflecting  boundaries  with  ^  0  which  adds  both  a  mirror  point  source  and  a  mirror 
line  source  to  satisfy  the  extrapolated  boundary  condition. “  In  the  case  of  =0,  their  result  is 

in  agreement  with  our  results.  Using  our  method,  corresponding  mirror  sources  can  also  be 
defined  for  the  cumulant  approximation. 

In  conclusion,  we  are  able  to  track  the  transport  of  photons  from  the  near  ballistic  to  the 
final  diffusive  stage  for  both  semi-infinite  and  slab  geometries  using  the  cTM.  Unlike  the  sDM 
and  cmDM,  we  can  obtain  an  accurate  picture  for  the  photon  migration  for  L  <  7/, ,  as  well  as  for 
L>11,.  For  t»ljc,  the  results  from  the  cTM  match  those  obtained  using  the  cmDM.  The 

effects  of  the  physical  boundaries  can  be  viewed  by  the  variations  in  the  amplitude  and  the  long- 
tail  decay  for  >  11, .  When  is  small,  the  cTM  provides  a  time-resolved  profile  that  can  be 

reasonably  compared  with  experimental  measurements  for  photon  pulse  propagation  in  clouds 
and  seawater. 

This  work  is  supported  in  part  by  Lockheed  Martin,  USAMRMC  (Award  DAMD17-98- 
1-8147)  DARPA,  NYSTAR,  and  NASA  IRA. 
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Figure  1.  The  moving  center  of  the  photon  density  function  and  the  diffusion  coefficients 
and  a  function  of  time  t.  For  t»  Ijc ,  stops  at  a  new 

position  located  exactly  /,  from  the  original  position. 

Figure  2.  Photon  density  for  an  infinite  medium  with  the  distance  between  the  source  and 
detector,  ,  is  given  by  (a)  =  2.01, ;  (b)  =  7.0/, ;  and  (c)  =  10.0/,  using 

/,  =100m,  n  =  1.33,  and =0. 

Figure  3.  For  semi-infinite  geometry,  locations  of  the  real  and  image  sources,  z^and  z,  , 
respectively,  for  (a)  the  CMDM  and  SDM  configurations  and  (b)  the  CTM 
configuration  if  the  real  source  is  located  at  z  =  0.  For  the  CTM  configuration,  the 

subscripts  i  and  /  associated  with  z,  and  z,  represent  the  initial  and  final  locations 
of  each  source.  An  extrapolated  boundary,  B,  is  located  at  a  distance  z  =  ccl,  from  the 
physical  boundary  of  the  semi-infinite  medium. 

Figure  4.  For  slab  geometry,  (a)  the  CMDM  and  SDM  configurations  and  (b)  the  CTM 
configuration  if  the  real  source  is  located  z  =  0 .  Extrapolated  boundaries,  Bi  and  Bi, 
are  located  at  distances  z  =  oil,  from  the  physical  boundaries  of  slab.  For  the  CTM 
configuration,  the  subscripts  i  and  /  associated  with  z^  and  z,  represent  the  initial 
and  final  locations  of  each  source. 

Figure  5.  For  the  SDM,  CMDM,  and  CTM,  the  photon  density  for  the  (a)  semi-infinite  and  (b) 
one-layer  uniform  slab  media  with  z^^^  =7.0/,  for  /,  =100m,  n  =  1.33,  andjU„  =  0. 

The  thickness  of  the  slab,  L ,  is  equal  to  z,d„  for  each  uniform  slab  media.  The 
results  for  Zj^^  =  2.0/,  are  the  same  as  those  for  the  infinite  medium  during  early  time. 
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List  of  Figures  cont’d 

Figure  6.  Using  the  CTM,  the  photon  density  for  the  infinite,  semi-infinite,  and  slab  media 
using  (a)  =2.0Z,  and  (b)  =1.01,  for  I,  =100m,  n  =  1.33,  and/i^  =0.  The 

thickness  of  the  slab,  L ,  is  equal  to  z^^^  for  each  uniform  slab  media. 
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Time-resolved  Fourier  optical  diffuse  tomography  is  a  novel  approach  for  imaging  of  objects  in  a  highly  scat¬ 
tering  turbid  medium  with  use  of  an  incident  (near)  plane  wave.  The  theory  of  the  propagation  of  spatial 
Fourier  components  of  the  scattered  wave  field  is  presented,  along  with  a  fast  algorithm  for  three-dimensional 
reconstruction  in  a  parallel  planar  geometry.  Examples  of  successful  reconstructions  of  simulated  hidden  ab¬ 
sorptive  or  scattering  objects  embedded  inside  a  human-tissue-like  semi-infinite  turbid  medium  are  provided. 
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1.  INTRODUCTION 

Research  on  the  use  of  near-infrared  diffusive  light  for 
biomedical  imaging  and  diagnosis  has  advanced  over  the 
past  decade  because  of  the  potential  of  the  technique  to  be 
a  safe,  noninvasive,  affordable,  and  superior  diagnostics 
In  the  search  for  a  methodology  that  provides 
fast  data  acquisition  and  reconstruction  to  perform  imag¬ 
ing  with  high  resolution  in  real  time,  a  variety  of  tech¬ 
niques  have  been  explored  including  the  use  of  time- 
resolved  picosecond  pulses,  continuous  waves,  and  diffuse 
photon-density  waves.  Most  methods  reconstruct  three- 
dimensional  (3D)  optical  property  maps  (OPMs)  by  matrix 
inversion,  by  iterative  techniques,  or  by  3D  rendering  of 
two-dimensional  {2D)  projection  images.^"®  The  degree 
of  difficulty  of  inverting  the  whole  3D  map  at  one  time  is 
usually  time  prohibitive  when  the  number  of  volume  ele¬ 
ments  involved  increases,  and  3D  rendering  of  two- 
dimensional  projection  images  requires  extra  depth  infor¬ 
mation  of  inhomogeneities  inside  turbid  media  to  behave 
well,  and  it  has  other  limitations.'^’® 

In  this  paper,  we  first  introduce  the  theory  of  propaga¬ 
tion  of  the  spatial  Fourier  component  of  the  scattered 
wave  field  inside  a  turbid  medium.  We  then  develop  a 
new  optical  diffuse  imaging  methodology  based  on  this 
theory,  using  the  two-dimensional  Fourier  transform  of 
photon  intensity  on  a  plane  to  detect  inhomogeneities  in  a 
highly  scattering  turbid  medium  when  the  medium  is  il¬ 
luminated  by  a  picosecond  (near-)  plane-wave  pulse.  In 
such  a  spatial  Fourier  space,  the  picture  of  photon  migra¬ 
tion  is  much  simplified  in  the  sense  that  different  spatial 
frequency  components  of  the  0PM  (2D  Fourier  transform 
on  the  xy  plane)  are  decoupled  from  one  another  and  de¬ 
pend  only  on  the  corresponding  spatial  frequency  compo¬ 
nent  of  the  photon  intensity  on  the  detector  plane.  On 
the  basis  of  this  observation,  we  obtain  a  super-fast  recon¬ 
struction  of  a  3D  0PM  by  matrix  inversion  of  each  spatial 
component  independently.  The  effect  of  noise  is  explic¬ 
itly  handled  by  controlling  the  set  of  spatial  frequency 
components  and  the  regularization  parameters  used  in 
the  matrix  inversion;  After  a  rigorous  account  of  the 


theory  and  a  brief  description  of  the  algorithm,  examples 
of  reconstruction,  by  using  backscattered  photons,  of  ab¬ 
sorptive  and  scattering  inhomogeneities  located  up  to  2 
cm  below  the  surface  of  a  human-tissue-like  semi-infinite 
turbid  medium  are  presented. 


2.  THEORY 

The  propagation  of  photon  density  (f>{r,t)  at  position  r 
and  time  t  in  a  highly  scattering  turbid  medium  can  be 
described  by  the  diffusion  equation 

d 

—  (^(r,  0  -  cV  •  Z)(r)V<^(r,  0  +  cpia{T^)^(T,t) 
dt 

=-S{r,t).  (1) 

The  absorption  coefficient  (per  unit  length),  and  the 
diffusion  coefficient  D  =  l/(3/xp,  where  /Xg  is  the  reduced 
scattering  coefficient,  may  depend  on  the  position  in  the 
medium;  c  is  the  speed  of  light  inside  the  medium,  and  S 
is  the  source  term  describing  the  density  of  photons  gen¬ 
erated  per  second. 

For  the  case  of  a  uniform  medium  and  an  incident 
source  S(r,t)  (S  =  0  when  ^  <  0),  the  incident  wave 
field  is  t)  =  /dV;od^'S(r',  ^ ')G(r,  r',  ^  -  V) 

where  G(r,  r',  t)  is  the  Green's  function  for  the  diffusion 
equation  in  a  uniform  turbid  medium.  When  some  weak 
inhomogeneities  (objects  such  as  tumors)  are  embedded  in 
the  medium,  we  vmte 

Ma,obj(r)  =  Ma  +  ^Aa(r), 

=  Ms  +  ^Ms(r),  (2) 

where  /x^  and  constant  absorption  and  reduced 

scattering  coefficients  of  the  otherwise  homogeneous  me¬ 
dium  and  /Xa,obj(r)  and  /Xg  obj(r)  are  the  absorption  and  re¬ 
duced  scattering  coefficients  of  the  embedded  inhomoge¬ 
neity  that  are  spatially  dependent.  Plugging  Eq.  (2)  into 
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Eq.  (1)  and  noting  the  diffusion  parameter  of  the  inhomo¬ 
geneity  Z)obj(r)  =  2)  +  SDir)  =  1/(3/^;)  -  Sfil(r)/ 

we  have 

d 

—  -  DcV^4>{r,t)  +  tiaC<t>{r,t) 

dt 


=  S(r,t)  +  cV  •  SD{r)V(i){r,t)  -  c  Wr)<A(i',  <).  (3) 

The  complete  right-hand  side  of  Eq.  (3)  now  acts  as  the 
source  term,  of  which  SiVyt)  contributes  to  the  unper¬ 
turbed  wave  field  <^o  ”  terms 

contribute  to  the  scattered  wave  field, 


t) 

=  <^(r,  t)  -  0o(r,  t) 


d^'G(r,  T\t  —  t*) 


X  [cVr^  .  8D(T^)V^.<t>{v\n 

-j  d^r'J‘dt’G(r,r’,t  -  t')S,i,(r')c,j>(r' ,t') 


+ 


Vr.G(r,  r',^  -  f') 


(4) 


after  integration  by  parts. 

To  first  order  in  the  variation  of  optical  absorption  and 
reduced  scattering  coefficients,  we  can  just  replace 
<f>{r\  in  Eq.  (4)  with  ,  i.e.,  the  total  wave  field  is  a 
superposition  of  the  incident  wave  field  <}>i  and  the  singly 
scattered  wave  field  .  This  is  the  well-known  Bom  ap¬ 
proximation. 

Now  consider  the  configuration  of  the  frequently  stud¬ 
ied  parallel  planar  geometry  (slab  or  semi-infinite)  with 
its  boundaries  at  2  =  0  and  z  —  d  {d  -  +00  for  semi¬ 
infinite  geometry).  The  exact  Green’s  function  is^® 


G(r,r\0 


1  /  ip-p'i" 


y.  G^{z,z',t),  (t>0),  (5) 

where  p  =  (*,3'),  p'  =  ix' ,y').  G^(z,z',t)  is  chosen 

according  to  the  boundary  condition  of  the  parallel  planar 
geometry  and  depends  only  on  time  t  and  the  z  coordi¬ 
nates  of  the  source  position  r  and  the  target  position  r'. 
Its  two-dimensional  Fourier  transform  on  p  is 


G{q,z,p',z',t)  =  J  d^pG(p,z,p',z',t)exp(-iq  ■  p) 

=  exp(— iq  •  p'  —.Dctq^ 

-  fi^ct)G2(z,z’,t) 

=  G(q,2,z',t')exp(-iq  •  p').  (6) 

For  simplicity,  we  restrict  ourselves  first  to  the  case  of 
a  pure  absorptive  perturbation  ^  0  ~  0) 


and  of  an  incident  pulse  S(r,t)  =  S(p)S(.z  -  Zs)S(.t). 
The  scattered  wave  field  on  a  plane  0  <  z  <  d  is  thus 

(f>siP,z,t) 

=  -  j  dV  J  d2p"J  dt'G(r,  r',t  -  t')SpLa(r') 

X  cS(p")G(r',p'',z,,t')  (7) 

firom  Eq.  (4)  after  ^  is  replaced  by  </>;.  Inside  Eq.  (7),  ex- 
p^md  the  source  term  S(p")  and  the  Green’s  functions 
G(r,  r',t  -  t')  and  G(r',  p",  z^ ,  f ')  into  integrals  of  their 
Fourier  components;  thus  we  find 

X  j  d^qG(q,z,z’,t  -  t') 

X  exp[jq  •  (p  -  p')] 

X  SiJi^(p',z’)c  I  d^q''S(q")exp(iq"  ■  p") 

X  j  d\G(q',z’,z,,t’) 

X  exp[iq'  •  (p'  -  p")] 

X  exp(iq  •  p)G(q, z,z',t  -  t') 

X  S(q")G(q',z',Zs,t') 

X  j  d^p’S/iia(p',z')exp[  -  ip' 

■  (q  -  q')]|  d2p"exp[ip''  •  (q"  -  q')] 

X  exp(8q  •  p)G(q,z,z',t  —  #') 

X  5p„(q  -  q',z'),§(q')G(q',z',Zs,t'), 

(8) 


where 

S(q)  =  S(q,zJ  =  j  d^pS(p,Zs)exp(-iq  ■  p), 

5/io(q.2)  =  j  d2p5yaa(p,z)exp(-jq  •  p) 

are  2D  Fourier  transforms  of  the  source  on  the  z  -  Zg 
plane  and  of  the  perturbation  of  the  absorption  coefficient 
over  the  2=2  plane,  respectively.  Finally,  we  recognize 
the  2D  Fourier  transform  of  the  scattered  wave  field 
<f>g{p,z,  t)  on  a  plane  2  for  the  case  of  a  pure  absorptive 
perturbation: 
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(j>{Ci,Z,t)  = 


dz'<5Aa(q~  q',^').S(q',^s) 


^s(q,2.0 


dz'  ^/ta(q.2')M'a(q.  0,Z,  t-,z') 


X  f  d/'G(q,z,2',<  -  /')G(q',z',z,,/'). 
Jo 

(9) 

In  a  similar  fashion,  for  the  case  of  a  pure  scattering 
perturbation  =  0  and  S/i'  0),  the  2D  Fourier 
transform  of  the  scattered  wave  field  is 


12irVs 

X  [‘dt’ 

Jo 


^2  d\'dz’S/i^(q  -  q’,z')S(q’,Zs) 


+ 


q  •  q'G(q,2:,2:',  t  -  t')G(q\  z' ,  ,  t') 

dG{q,z,z\t  -  V)  dG{q:,z\z,,t^) 


dz’ 


dz' 


(10) 


The  general  Fourier  scattered  wave  field  is  the  sum  of 
Eq.  (9)  and  Eq.  (10).  Denoting  the  convolutions 

^  \  dt'G{q,z,z\t 
Jo 

-  t')G{q\z\Zs,t'), 

n  dG{qyZ,z' ,  t  -  t') 
Ws(ci,q'yZ,t;z')  =  d^' — 

Jo 


fo  dz' 

dG{q'  ,z'  ,Zs,t') 


X 


dz' 


(11) 


which  are  the  weight  functions  involved  in  the  propaga¬ 
tion  of  spatial  Fourier  components  of  the  scattered  wave 
field,  we  have 


=  J  d^q'dz' 

,  c 

X  S{q! ,Zs)Wa{ii.^! +  ^27rVs^ 

X  J  d^q' dz' dji'^iq  -  q', 2')»S(q', Zg) 

X  [q  •  q'M;a(q»  q! ,z,t\z')  +  z^s(q>  ^z,t\z')']. 

(12) 

For  the  simple  case  in  which  the  incident  wave  is  a 
plane-wave  pulse  (see  Appendix  A  for  justification),  i.e., 
S(r,t)  =  Sd{z  -  Zs)S{t)  where  S  is  a  constant,  such 
thatS(qyZs)  =  47T2S5(q),  Eq.  (12)  simplifies  to 


sK(^,z') 

3m;" 


u;5(q,  0,z,t;z') 


(13) 


The  most  salient  feature  of  the  above  result  [Eq.  (13)]  is 
that  different  spatial  frequency  components  of  Sfia  and 
Sfi's  are  decoupled  from  one  another,  and  the  q  component 
of  the  optical  parameters  depends  only  on  the  correspond¬ 
ing  spatial  frequency  component  of  the  scattered  wave 
field  ^siq.Zyt).  Thus  the  dimension  of  the  inverse  prob¬ 
lem  to  be  solved  below  is  greatly  reduced,  as  is  the  com¬ 
putation  time. 

If  we  approximate  the  integration  over  2'  by  a  summa¬ 
tion  and  fix  2:  =  2:^  at  the  detection  plane  (omitting  z  here¬ 
after),  the  Fourier  scattered  wave  field  on  the  detection 
plane  is 


^,(q,  t)  =  Sc^z^ 
>=i 


-^Aa(q.  27)«’a(q. 


+ 


w^(q,(i,t-,Zj)  , 


(14) 


where  Az  is  the  discretized  step  size,  is  the  total  num¬ 
ber  of  slices  (layers)  in  the  z  direction  between  the  source 
plane  and  the  detection  plane,  and  Zj  is  the  z  coordinate  of 
the  central  position  of  layer  7. 

If  we  set  q  =  0  in  Eq.  (14), 


^,(0,/)  =  ScAz^ 

/=i 


-Sjl^{0,Zj)w^(0, 0,t;zj) 


+ 


— — u;s(0,0,f;zy)|, 


(15) 


the  zero  spatial  frequency  components  5/i-a(0,Zy)  and 
Sfig{0,Zj)  can  be  readily  solved  without  the  need  for  a 
complete  reconstruction.  Owing  to  the  nature  of  Fourier 
transform,  they  just  provide  the  profile  of  the  amount  of 
total  perturbation  of  absorption  and  reduced  scattering 
coefficients  per  slice,  i.e.,  the  depth  profile  of  the  inhomo¬ 
geneities. 

The  whole  3D  map  of  absorption  and  reduced  scatter¬ 
ing  coefficients  is  thus  constructed  through  an  inverse 
Fourier  transform  from  all  the  q  components  of  Sfia  and 
Sfll  at  different  depths,  each  of  which  is  solved  indepen¬ 
dently  from  a  series  of  time-resolved  scattered  wave  field 
by  Eq.  (14). 

A  schematic  diagram  of  the  procedure  of  image  recon¬ 
struction  is  shown  in  Fig.  1.  The  maximum  spatial  fre¬ 
quency  (cutoff  frequency)  of  the  components  used  in  the 
inversion  is  determined  through  a  signal-to-noise-ratio 
analysis  in  which  the  Fourier  components  whose  magni¬ 
tudes  fall  below  a  threshold  (comparable  to  the  noise 
level)  are  discarded.  The  regularization  parameter  in 
the  matrix  inversion  is  obtained  by  the  robust  L-curve 
method.  The  L-comer  finder,  which  locates  the  comer 
by  maximum  curvature, is  implemented  and  is  used  to 
obtain  the  regularization  parameter.  Neither  visual  es- 
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Fig.  1.  Schematic  diagram  of  image  reconstruction. 


Fig.  2.  Geometry  for  time-resolved  Fourier  optical  diffuse  to¬ 
mography  with  use  of  backscattered  photons.  The  source  is  a  pi¬ 
cosecond  (near-)  plane-wave  pulse  and  a  series  of  snapshots  of  a 
10  X  10  cm^  area  on  the  surface  are  computed  as  the  input  to  im¬ 
age  reconstruction.  Absorptive  objects  A  (—2.5,  —1.875,  -0.75) 
cm,  B  (-1.25,  -0.31,  -0.75)  cm,  C  (0.94,  1.56,  -1.95)  cm,  and  D 
(0,94,  -0.625,  -1.95)  cm  or  scattering  objects  E  (-2.5,  -1.875, 
-0.75)  cm,  F  (-1.25,  -0.31,  -0.75)  cm,  G  (0.94, 1.56,  -1.35)  cm, 
and  H  (0.94,  —0.625,  —1.35)  cm  are  used  in  the  simulation. 


timate  nor  prior  information  is  required  for  this  proce¬ 
dure.  L  curves  are  different  for  each  spatial  frequency  q. 
The  regularization  parameter  is  determined  from  the  re¬ 
construction  of  depth  profile  (v^here  an  inversion  for  q 
=  0  is  performed).  The  same  value  is  then  used  in  the 
full  3D  reconstruction  (layer  reconstructions,  where  in¬ 
version  includes  q  0). 

Both  transmission  and  backscattering  image  recon¬ 
struction  configurations  can  easily  be  made  by  using  Eqs. 
(13)  and  (14). 


3.  SIMULATION 

For  demonstration  purposes,  consider  a  semi-infinite  tur¬ 
bid  medium  (2:  <  0)  with  its  surface  at  2;  =  0  (Fig.  2), 
whose  absorption  coefficient  =  0.0033  mm“^  and  re¬ 
duced  scattering  coefficient  =  1.0  mm“^. 

A.  Absorptive  Inhomogeneity 

Four  absorbing  objects  A,  B,  C,  and  D,  each  6.25  mm 
X  6.25  mm  X  3  mm  and  with  absorption  coefficient 
/^c,obj  “  0.02  mm“^  and  reduced  scattering  coefficient 
equal  to  that  of  the  background,  are  placed  at  depth  7.5, 
7.5,  19.5,  and  19.5  mm  below  the  surface,  and  their  xy  co¬ 
ordinates  are  (-25,  —18.75),  (-12.5,  -3.1),  (9.4,  15.6), 
and  (9.4,  6.25)  mm,  respectively.  The  medium  is  illumi¬ 
nated  by  an  incident  pulse  of  a  Gaussian  shape  of 
exp(-pV2cr^)  with  a  —  50  mm  inside  an  aperture  of  ra¬ 
dius  50  mm,  propagating  along  the  negative  z  axis  at  time 
t  =  0. 

These  parameters  are  potentially  applicable  to  optical 
mammography  of  the  compressed-breast-toward-chest 
setup  with  use  of  backscattered  photons.  A  series  of 
simulated  measurements  (total  =  15  snapshots  from 
300  to  2400  ps)  of  an  area  100  mm  X  100  mm  on  the  sur¬ 
face  plane  2:  =  0  are  generated  by  using  a  direct  calcula¬ 
tion  for  the  Gaussian  pulse  in  r  space.  The  simulated 
data  are  used  as  input  for  inversion  after  adding  a  1%, 
5%,  or  10%  Gaussian  noise. 

In  the  reconstruction  part,  the  near-surface  region 
of  the  turbid  medium  of  depth  up  to  3  cm  is  sliced  into 
Nz  =  10  layers,  i.e,,  Az  =  0.3  cm,  and  objects  A  and  B  are 
then  located  on  layer  3,  and  C  and  D  are  located  on  layer 
7.  The  detection  plane  of  an  area  of  10  X  10  cm^  is  di¬ 
vided  uniformly  into  a  NJSly  =  32  X  32  grid.  Objects  A, 

B,  C,  and  D  all  take  2X2  elements  by  this  grid.  The 
results  of  reconstruction  are  shown  below. 


(a)  (b)  (c) 

Fig.  3.  Absorption  depth  profile  for  (a)  with  1%  noise,  (b)  5%  npise,  and  (c)  10%  noise. 
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1 .  Depth  Profile 

The  absorption  depth  profile,  i.e.,  the  total  absorption  per¬ 
turbation  per  layer  /d^p^/ta(p,  2)  versus  depth  z  is  shown 
in  Fig.  3  with  different  noise  levels  for  cases  (a)  1%  noise, 
(b)  5%  noise,  and  (c)  10%  noise.  In  Fig.  3(a)  there  is  one 
peak  at  depth  2  =  0.75  cm  (layer  3)  where  objects  A  and  B 
are  embedded,  and  another  peak  at  2  =  1.95  cm  (layer  7) 
where  objects  C  and  D  are  embedded.  The  width  of  the 
first  peak  at  half-height  is  0.34  cm,  approximately  the 
thickness  of  one  layer  (0.3  cm),  which  means  that  the 
depth  of  objects  A  and  B  is  resolved  very  well.  The  sec¬ 
ond  peak  of  objects  C  and  D  spans  two  and  a  half  layers 
with  its  width  of  peak  at  half-height  0.74  cm,  but  its  peak 
position  is  still  correct. 

When  the  level  of  noise  increases,  the  peak  values  of 
both  peaks  decrease,  and  the  half-widths  increase.  The 
effect  on  the  second  peak  at  2  =  1.95  cm  is  much  more 
significant  than  that  on  the  first  one  at  2  =  0.75  cm. 

2.  Layer  Reconstruction 

The  full  3D  0PM  is  reconstructed.  The  reconstructed  ab¬ 
sorption  coefficients  of  the  layers  at  the  two  peak  posi¬ 
tions  are  shown  in  Figs.  4-6  for  the  three  noise  levels. 


Figure  4  shows  that  objects  A  and  B  are  clearly  resolved 
as  two  objects  centered  at  their  original  positions  with 
negligible  expansion;  and  objects  C  and  D  at  depth  2 
=  1.95  cm  are  also  detected  at  the  correct  central  posi¬ 
tions,  but  the  resolved  images  are  expanded  on  the  xy 
plane.  With  an  increase  in  noise  level,  the  shape  of  ob¬ 
jects  A  and  B  blurs  from  Figs.  4(a)  to  5(a)  and  6(a),  and 
the  blur  is  even  worse  for  objects  C  and  D  under  the  same 
condition  [from  Figs.  4(b)  to  5(b)  and  6(b)]. 

At  noise  level  of  1%,  the  reconstructed  absorption  pa¬ 
rameter  for  objects  A  and  B  is  0.0071  mm”^  approxi¬ 
mately  36%  of  the  original  value  0.02  mm“^  of  the  absorp¬ 
tive  inhomogeneity.  In  other  words,  the  object  appears 
larger  in  space  with  a  weakened  absorption  parameter. 
As  the  noise  level  increases,  the  effect  is  accentuated  with 
a  further  reduction  in  the  resolved  absorption  parameter. 

B.  Scattering  Inhomogeneity 

For  another  example,  four  scattering  objects  E,  F,  G,  and 
H,  each  6.25  mm  X  6.25  mm  X  3  mm  and  with  reduced 
scattering  coefficient  obj  =  ^.5  mm“^  and  absorption  co¬ 
efficient  equal  to  that  of  the  background,  are  placed  at 
depth  7.5,  7.5,  13.5,  and  13.5  mm  below  the  surface,  and 


0.001 


0 


(a)  (b) 

Fig.  4.  Layer  reconstruction  at  a  noise  level  of  1%:  (a)  resolved  objects  A  (left)  and  B  (right)  at  2  =  0.75  cm  (layer  3);  (b)  resolved 

objects  C  (upper)  and  D  (lower)  at  2  =  1.95  cm  (layer  7).  The  darkness  of  the  pixel  represents  the  resolved  absorption  coefficient  in  units 
of  inverse  millimeters. 


X(cm)  X(on) 

(a)  (b) 

Fig.  5.  Layer  reconstruction  at  a  noise  level  of  5%:  (a)  resolved  objects  A  (left)  and  B  (right)  at  2  =  0.75  cm  (layer  3);  (b)  resolved 

objects  C  (upper)  and  D  (lower)  at  2  =  1.95  cm  (layer  7).  The  darkness  of  the  pixel  represents  the  resolved  absorption  coefficient  in  units 
of  inverse  millimeters. 
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(a)  (b) 

Fig.  6.  Layer  reconstruction  at  a  noise  level  of  10%:  (a)  resolved  objects  A  (left)  and  B  (right)  at  z  =  0.75  cm  (layer  3);  (b)  resolved 
objects  C  (upper)  and  D  (lower)  at  z  =  1.95  cm  (layer  7).  The  darkness  of  the  pixel  represents  the  resolved  absorption  coefficient  in  units 
of  inverse  millimeters. 


(a)  (b)  (c) 

Fig.  7.  Scattering  depth  profiles  for  (a)  with  1%  noise,  (b)  5%  noise,  and  (c)  10%  noise. 


(a)  (b) 

Fig.  8.  Layer  reconstruction  at  a  noise  level  of  1%:  (a)  resolved  objects  E  (left)  and  F  (right)  at  z  =  0.75  cm  (layer  3);  (b)  resolved  objects 
G  (upper)  and  H  (lower)  at  z  =  1.35  cm  (layer  5).  The  darkness  of  the  pixel  represents  the  resolved  reduced  scattering  coefficient  in 
units  of  inverse  millimeters. 


their  xy  coordinates  are  (-25,  —18.75),  (—12.5,  —3.1), 
(9.4,  15.6),  and  (9.4,  6.25)  mm,  respectively.  Objects  E 
and  F  are  now  located  on  layer  3,  and  G  and  H  are  located 
on  layer  5.  The  same  source  and  inversion  procedure 
used  in  the  previous  example  are  used  here.  The  results 
of  reconstruction  are  shown  below. 


1.  Depth  Profile 

The  scattering  depth  profile  is  shown  in  Fig.  7  with  differ¬ 
ent  noise  levels  for  cases  (a)  1%  noise,  (b)  5%  noise,  and  (c) 
10%  noise.  Two  peaks  are  correctly  revealed  with  the 
first  at  depth  z  =  0.75  cm  (layer  3)  and  another  at  z 
=  1.35  cm  (layer  5),  where  objects  E  and  F,  G  and  H  are 
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embedded,  respectively.  In  Fig.  7(a),  the  width  of  the 
first  peak  at  half-height  is  0.3  cm  and  that  of  the  second 
peak  is  0.42  cm.  This  means  that  the  depth  of  these  ob¬ 
jects  is  resolved  quite  well. 

With  an  increase  in  noise  level,  we  observe  the  same 
behavior  of  decreasing  quality  of  depth  resolution  as  in 
the  absorptive  case. 

2.  Layer  Reconstruction 

The  reconstructed  scattering  coefficients  of  the  layers  at 
the  two  peak  positions  (layer  3  and  layer  5)  are  shown  in 
Figs.  8-10  for  the  three  noise  levels. 

We  observe  a  result  similar  to  that  for  the  absorptive 
case.  Objects  E  and  F  are  better  resolved  than  objects  G 
and  H,  which  are  deeper  into  the  turbid  medium,  and  the 
noise  has  a  more  adverse  effect  on  objects  G  and  H  than 
on  objects  E  and  F.  The  reconstructed  reduced  scattering 
coefficient  of  objects  E  and  F  is  0.27  mm"\  approximately 
54%  of  the  original  value  of  the  embedded  scattering  in¬ 
homogeneity,  at  a  noise  level  of  1%. 


4.  DISCUSSION 

A  fast  time-resolved  Fourier  optical  diffuse  tomography 
based  on  decoupled  propagation  of  spatial  Fourier  compo¬ 
nents  of  the  scattered  wave  field  when  the  medium  is  il¬ 
luminated  by  a  plane  wave  is  presented.  For  a  wave  not 
strictly  plane  but  whose  zero-frequency  component  domi¬ 
nates,  this  approximation  is  still  valid  as  long  as  the  ra¬ 
dial  dimensions  of  the  volume  where  inhomogeneities  ex¬ 
ist  are  much  smaller  than  the  effective  width  of  the 
Gaussian  beam  (see  Appendix  A). 

The  image-reconstruction  method  provided  is  efficient. 
The  optical  parameters  atN^NyN^  different  voxels  are  re¬ 
constructed  from  a  set  ofiV^  measurements  hyNj  times  of 
inversions  of  X  matrices,  where  NI  <  N^Ny  is  the 
total  number  of  Fourier  components  with  the  noisy  ones 
discarded.  Our  procedure  is  much  more  efficient  than  a 
direct  reconstruction,  where  an  inversion  of  N^NyN^ 
X  NJ^yN^  matrix  is  involved.  The  speedup  is  approxi¬ 
mately  OiN^Ny)  times  faster.  The  time  this  algorithm 
takes  to  perform  a  complete  3D  reconstruction  in  the 


(a) 

Fig.  9.  Layer  reconstruction  at  a  noise  level  of  5%: 
G  (upper)  and  H  (lower)  at  z  =  1.35  cm  (layer  5). 
units  of  inverse  millimeters. 


(b) 

(a)  resolved  objects  E  (left)  and  F  (right)  atz  =  0.75  cm  (layer  3);  (b)  resolved  objects 
The  darkness  of  the  pixel  represents  the  resolved  reduced  scattering  coefficient  in 
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(a) 

Fig.  10.  Layer  reconstruction  at  a  noise  level  of  10%: 
objects  G  (upper)  and  H  (lower)  at  2  =  1.35  cm  (layer  5). 
in  units  of  inverse  millimeters. 


(b) 

(a)  resolved  objects  E  (left)  and  F  (right)  at  z  =  0.75  cm  (layer  3);  (b)  resolved 
The  darkness  of  the  pixel  represents  the  resolved  reduced  scattering  coefficient 
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above  examples  (of  32  X  32  X  10  volume  elements)  is  less 
than  half  a  minute  with  use  of  the  scripting  language  Py¬ 
thon  on  one  180-Mhz  CPU  of  an  Origin  200  computer 
from  Silicon  Graphic  Inc.  This  algorithm  scales  only  lin¬ 
early  with  the  number  of  elements  in  the  xy  grid,  so  it  can 
be  used  to  handle  much  larger  data  sets  in  real  time  with 
little  difficulty. 

This  approach  does  not  limit  the  number  or  the  thick¬ 
ness  of  the  inhomogeneities.  It  allows  multiple  inhomo¬ 
geneities,  and  one  inhomogeneity  may  span  several  lay¬ 
ers. 

With  little  effort,  a  depth  profile  (the  sum  of  the  pertur¬ 
bation  of  the  optical  parameter  versus  depth)  of  the  inho¬ 
mogeneities  inside  a  highly  scattering  turbid  medium  can 
be  obtained.  This  information  itself  may  be  very  useful 
in  some  cases.  When  the  inhomogeneity  is  found  to  exist 
only  in  one  layer  from  the  depth  profile,  the  summation  in 
Eq.  (14)  no  longer  exists.  A  direct  inverse  Fourier  trans¬ 
form  can  thus  be  used  to  resolve  the  inhomogeneity  when 
it  is  a  solely  absorptive  or  scattering  perturbation. 


APPENDIX  A 

Equation  (12)  is  the  exact  formula  for  calculating  the 
scattered  wave  field.  For  a  pulse  Sir,  t)  =  Sip) 
Siz  -  Zs)Sit)  with  Gaussian  shape  Sip) 

~  Soexp(-pV2c7^),  we  have  Siq)  =  27ro-^So 
X  expi-a^qV2)  and  the  first  term  of  Eq.  (12), 

cr^SnC  f 

A  =  — - -  dz' S/Maifl  -  q',2') 

277  J 

X  ex-p{-a^q'V2)Wa{q.,  q',z,t;z') 
a^SoC  f 

=  -  d‘^q'dz'd^p'dn,ip',z') 

277  J 

X  exp[-i(q-  q')  •  p'  -  a^q’V2] 

X  Wa{q,q',z,t;z') 

Cr^SnC  f 

=  — r -  d^p'dz'Sfiaip',z')exp{-iq  ■  p') 

277  J 

X  j  dt'Giq,z,z\t  ~  t') 

Jo 

X  GJz',Zs,t')expi-/iaC^') 

X  J  d^q'  exp(-o-^q'^/2  -  Bet'  +  iq'  •  p') 

(16) 

after  plugging  in  Eq.  (11)  and  Eq.  (6). 

The  last  integral  of  Eq.  (16)  can  be  performed  exactly 
and  turns  out  to  be  7rB~^  exp(-4p'^/i?^),  where  the  effec- 
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tive  width  ^  a ^12  +  Dct'  >  a^l2.  When  the  inho¬ 
mogeneities  exist  inside  a  region  of  radial  dimension  L 
around  the  origin  of  the  xy  coordinate  system  that  satis¬ 
fies  L  <  R,  we  can  approximate  exp(-4p'^/i2^)  by  1, 
which  is  equivalent  to  letting  q'  — >  0,  the  case  of  an  inci¬ 
dent  plane  wave.  The  error  made  by  such  an  approxima¬ 
tion  is  of  second  order  in  L/R . 

The  same  analysis  can  be  applied  to  the  second  term  of 
Eq.  (12). 
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